Introduction
What is a partially ordered set?
A partially ordered set (or poset) is a collection of things where we have a nice way to say that some things are smaller or bigger than other things, but we don't necessarily have a way to compare every pair of elements, and this way of comparing things makes sense. The most basic example would be to take a subset of the integers, and we have an intuitive notion that when we take a pair of numbers, one is going to be smaller than (or equal to) the other. That would be a totally ordered set. But there are other examples!
One non-trivial but still elementary example would be to take all the divisors of a positive integer, and compare them by seeing if one evenly divides into a another. For example, the divisors of 12 are 1,2,3,4,6,12. We can see that we have relations like 1 is a divisor of 3, and 2 is a divisor of 4. But we can't compare everything. 3 isn't a divisor of 4, and 4 isn't a divisor of 3. If this ever comes up again, we'll call it the divisor poset.
Another non-trivial but still elementary example would to take the set of subsets of a given set, and see if one subset contains another. For example, if we had the set {a,b,c}, then we would have the subsets $\emptyset$,{a},{b},{c},{a,b},{a,c},{b,c},{a,b,c}. Again, we can compare some things ({a} is a subset of {a,c}) but not other things ({b} is not a subset of {a,c}, nor vice versa). If this every comes up again, we'll call it the boolean poset.
*Sidebar: If you want to search for these things, they're usually called divisor lattice and boolean lattice. A lattice is just a poset with some additional special properties.
When we want to talk about partially ordered sets without a particular example in mind, we'll typically call the set $P$ and use the relation $\leq_P$. The basic things we need to have a nice/sane notion of "smaller/bigger" are that we always have $x\leq_P x$ (reflexivity), if $x\leq_P y$ and $y\leq_P x$ then $x=y$ (anti-symmetry) and if both $x\leq_P y$ and $y\leq_P z$, then we must have $x\leq_P z$. These are some basic things we would expect a notion of "smaller/bigger" to satisfy, and keeps us from making loops. For brevity, unless we really need to distinguish between a poset relation $leq_P$ and the usual $\leq$ relation, we will drop the subscript $P$.
A convenient way to draw/visualize a poset is to think of the elements as nodes in a graph, and we have an edge going from $z$ down to $x$ if $x\leq z$, and there are no other elements $y$ with $x\leq y\leq z$. This is called a covering relation, and we sometimes denote it $x\lessdot y$. We call resulting graph the Hasse diagram of $P$, and often use it interchangeably with $P$.
This is very close to the notion of an acyclic digraph. A key technical difference is that if we have a<-b and b <- c in a digraph, we may or may not have the edge a<-c . For a poset, that additional edge would never appear in the Hasse diagram, as the transitivity property means we know that relation always holds. A more stylistic difference is that we don't use directed arrows in a Hasse diagram, it is just always implied that things lower in the diagram/graphic are smaller in the poset.
EXAMPLES
What is a linear extension?
A linear extension of a partially ordered set is simply a way of making a bigger set of relations such that any relations that previously held are still true. We are extending our relation to something where everything can be arranged in a linear fashion.
Take the divisors of 12, as an example. Then we can see that our normal way of comparing two integers (in terms of size) is a linear extension. If we say $1\leq 2\leq 3\leq 4 \leq 6\leq 12$ , then it's true that anytime we had $x | y$, then it will be true that $x\leq y$. But there are many other ways we could order these numbers are still have that be true. We could order them 1,3,2,4,6,12. Or 1,2,4,3,6,12. As an exercise, see if you can list out all 5 linear extensions.
This is sometimes the more traditional way to think of a linear extension, as it's a simple array of your poset elements. But I prefer (and think it's more instructive) to think of is as a way of numbering nodes in the Hasse diagram so that nodes with smaller numbers always appear lower. It's also closer to the (essentially equivalent) notion of a topological sort for an acyclic directed graph.
EXAMPLES
Say we wanted to know how many linear extensions there are. One naive way to count the number of linear extensions of a partially ordered set would be to just generate all of them and count how many you get. There are nice ways to generate them, one at a time, making sure you don't miss any. But the number of linear extensions of a poset can be extremely large even for posets with a small number of elements. So is there a way that we can count without enumerating?
What are order ideals?
An order ideal of a poset is a subset of elements that is "closed under going down". That is, $I$ being an order ideal means if $y\in I$ and $x\leq_P y$, then $x\in I$.
EXAMPLE So if we had a poset $P$ with a<b, a<c, b<d, c<d (your Hasse diagram should look like a diamond), the order ideals would be empty set, {a}, {a,c}, {a,b}, {a,b,c}, and {a,b,c,d}. Something like {b,c} wouldn't be an order ideal, because $b\in I$, and $a\leq_p b$, but $a\not\in I$. It's pretty easy to see that we'll always at least have the empty set and the whole set as order ideals.
In a bit of a case of inception, these order ideals themselves are a partially ordered set(!), using the relation of set containment. We call this the lattice of order ideals, and denote it $J(P)$. EXAMPLE
One notable example is if we have a poset with no relations between distinct elements, just $x\leq_P x$, and this is called the 'antichain poset'. Then every subset is an order ideal, and $J(P)$ ends up being a boolean poset.
As a side note, there's a pretty cool theorem (Birkhoff's theorem of finite distributive lattices) that says if the relations of your poset act enough like set containment, then you can always construct as associated poset whose lattice of order ideals is isomorphic to your original poset.
Why do we care about order ideals?
The short answer is that "linear extensions in $P$ are maximal chains in $J(P)$, and maximal chains are easy to count without generating all of them".
First, we see that linear extensions in $P$ are in fact equivalent to maximal chains in $J(P)$. Start off with some linear extension, think of it as a labelling of the elements with $1,\ldots,|P|$. Look at the set of elements with labels <1. It's empty (dumb corner case we need to start with). Then look at the set of elements with labels less than 2 (just the element labeled 1). Then less than 3. And so on, until we have the whole poset. At every step, we must have an order ideal. Say we're looking at everything with a label less than $k$. If $y\in P$, then $y$ has a label less than $k$. If $x\leq_P y$, then $x$ has a label smaller than $y$ So $x$ must have a label smaller than $k$, and $x\in I$. Going in reverse, we can see that everytime we make a step up in a maximal chain of $J(P)$, we're adding a new element, and so we just number the elements we add in the order we add them.
In terms of maximal chains being easy to count, it's a standard trick in graph-theory-esque problems of breaking up a path from a to b into paths from a to things one step away from b, and then the last step. In this case, the number of saturated chains from a minimal element of our poset to a given element $x$ is the sum of the number of saturated chains from a minimal element of our poset to each element that is covered by $x$. We're a little fortunate because we're looking at maximal chains in $J(P)$, which always has a unique minimal element.
EXAMPLE Say we have the boolean lattice on $\{a,b,c\}$. The empty set is the only minimal element, only one path to itself. For each of the singletons ({a},{b},{c})), there is again only one saturated chain from the bottom to that element. For each of the doubletons ({a,b},{a,c},{b,c}), there are two paths from the bottom to each element. And then for the maximal element, there are 6 paths. They key point is that to figure out the maximal element had 6 chains, we didn't really need to go back and look at all 6 chains. We could have just known that it covers $\{a,b\},\{b,c\},\{a,c\}$, previously recorded that each of those things had 2 chains going up to them, and adding two together three times to get 6.
Sidenote: In this particular example with the boolean lattice, it's easy to see that maximal chains in the boolean poset (or equivalently, linear extensions of a set where the elements have no non-trivial relations) are exactly in bijection with permutations. So in this case, we could have gotten $3!$ without any work. Typically, we are not so fortunate.
Why did I want to code this?
Historically, computational tools for mathematics software has had some downfalls. Either you hoped that Matlab/Maple/Mathematica would include a feature you wanted (and had no particular way to request or implement yourself), and have it be bug-free, and have it not change behavior across releases. Or, you could have somebody write custom packages on top of these software systems, which are not always necessarily well-documented, nor necessarily maintained across version updates, and it's difficult to add any improvements/modifications that may be beneficial without the help of the original author. And this is all assuming you were affiliated with a university that had an affiliation with the correct software provider.
In terms of specifics for why I wanted to code this, the algorithm originally comes from a Maple package for posets written by John Stembridge ( http://www.math.lsa.umich.edu/~jrs/maple.html#posets ). To my knowledge, the code still works if you have access to Maple, is extremely efficient, and is actually very well documented from a user perspective. But it does require you to have access to Maple, the code is not necessarily documented from an algorithmic perspective, and part of the reason for it's speed is the limitation that you only work on partially ordered sets whose elements are positive integers. By including this code as part of the SageMath project, it is accessible to everybody as part of an open-source software system, it will be maintained over time, there is better documentation for the algorithm, and a user doesn't have to do any work trying to make a naturally labelled poset on positive integers that is isomorphic to whatever they're interested in.
Naturally Labelled Posets
We say that a partially ordered set of $n$ elements is naturally labelled if the elements are $1,\ldots, n$, and if we label an element with itself, we get a linear extension. This can be a bit confusing from the approach of thinking about linear extensions as a labelling of the Hasse diagram. For example, we could have (1)<(2), (1)<(3), (2)<(4), and (3) < (4) as a naturally labelled poset, because if I gave (1) the label 1, (2) the label 2, (3) the label 3 and (4) the label 4, it would form a linear extension. But labelling (1) as 1, (2) as 3, (3) as 2, and (4) as 4 would also work. I've included parentheses to distinguish between poset elements and labels here, but there's not much of an equivalent when coding with lists besides making sure you're keeping track of what's an element and what's a label. Notably, this confusion doesn't arise thinking of linear extensions as an ordered list. For the previous example, we would have $[1,2,3,4]$ and $[1,3,2,4]$, and the labels are hidden as indices in the array.
A major benefit of restricting to working with naturally labelled posets is that you can do things like use a dictionary where the value for a key $x$ is the set of elements that cover $x$ (ie, a table of up-edges in the Hasse diagram). If your poset element $x$ is something like a subset in the boolean poset, it's not hashable, no dictionary.
The other benefit to working with naturally labelled posets is that if we traverse the elements in order (1,2,...), we are assured that we're going to be working our way from the bottom up. That means every time we reach an element, we know that we have already seen every element that it covers.
The big key to making this code work efficiently is that if we try to construct $J(P)$ directly from the definition, we are going to be getting elements that are subsets, and not positive integers. Since we really just need to count chains in $J(P)$, the magic is going to come by constructing a poset that's isomorphic to $J(P)$, but a naturally labelled one. As a lovely added detail I hadn't noticed until writing this, the element $k$ in our naturally labelled version of $J(P)$ will correspond to the $k^{th}$ order ideal of $P$ in lex-order.
Data Structure
We will represent our poset as a dictionary of up-edges in the Hasse diagram. That is to say, $P[x] = \{y\in P | x\lessdot y\}$. As an example, if we had the diamond poset $0\leq 1$, $0\leq 2$, $1\leq 3$, $2\leq 3$, then we would input $P$ as $\{0:[1,2], 1: [3], 2: [3], 3: []\}$. Note the zero-indexing for the input.
In Stembridge's posets package, one could also represent the poset as a list of the covering relations given as an ordered pair ($P = [[1,2],[1,3],[2,4],[3,4]]$), or a list of all the relations ($P = [[1,2],[1,3],[1,4],[2,4],[3,4]]$) (this implementation was not zero-indexed).
In SageMath, one could additionally give the input as a set and a function/lambda operator that returns a boolean as to whether $x\leq_P y$.
In any case, the dictionary/table of up-edges will be eventually calculated, so for the sake of having this be a self-contained notebook, we'll take that as our starting point.
P={}
P[0]=[1,2]
P[1]=[3]
P[2]=[3]
P[3]=[]
P
Our first step will be to convert our input of a dictionary of up-edges into a dictionary of up-relations. That is, instead of having $P[x] = \{y| x\lessdot y\}$, we will be working with $up[x] = \{y| x\leq y\}$. For some additional terminology, similar to how we had an order ideal being a subset closed under "going down", this will be an order filter (subset closed under going up). When we take all things smaller/bigger than a single element, that's called the principal order ideal/filter (resp.).
We set $n$ to be the number of elements in the poset.
We take advantage of the fact that the order filter of $x$ is the union of the order filters for elements that cover $x$, and reverse indices from top-down to ensure we've already computed the full order filter for every element covering $x$ before we reach $x$.
n = len(P.keys())
up = {key: value[:] for key, value in P.items()}
for i in range(n):
up[n - 1 - i] = sorted(set(up[n - 1 - i] + [item for x in up[n - 1 - i] for item in up[x]]))
Next, we initialize an empty dictionary for the up-edges in $J(P)$. We will use $m$ as a counter for how many elements we have currently added to $J(P)$. It seems that in my implementation, I did not have $J$ set to be zero-indexed. It was important for the input to be zero indexed, to match the fact that SageMath spits out a zero-indexed dictionary of edges from the Hasse diagrm for a poset. For this application, it is less important for the output to be zero-indexed to match Python convention, though I may correct this later.
We will perform a loop $n$ times, where after $k$ times running the loop, Jup will be the dictionary of up-edges for the poset of order ideals of $P$ restricted to the elements $0,\ldots, k-1$.
An additional piece of information we need is an auxilliary array $loc$. After we have finished $k$ loops, $loc[j]$ will tell us which element in the partially constructed poset of order ideals corresponds to the principal order ideal for $j$ (restricted to the elements $0,\ldots k-1$).
EXAMPLE
Say we have ran the loop twice. At this point, in $Jup$, we would have 1 (corresponding to the empty order ideal), 2 (corresponding to the order ideal $\{0\}$), and 3 (corresponding to the order ideal {0,1}). Then the array loc would be $[2,3,2,3]$.
Jup = {1: []}
m = 1
loc = [1] * n
Now, we consider what happens when we introduce a new element $x$ to the poset, and how it affects $Jup$. Recall that $loc[x]$ is the principal order ideal for x in the poset restricted to $0,\ldots, x-1$. Any new order ideal created must contain $x$, and so must contain $loc[x]$, and in general will be something above $loc[x]$ in $Jup$ but with $x$ added.
So we construct the set $K$ corresponding to the elements of $Jup$ that are above $loc[x]$. We initialize $K$ with the singleton list $loc[x]$, and then we append the list of all upper covers of $loc[x]$, and then all the upper covers of everything in that set, and so on, until we stop adding anything. Then we flatten this into a single list, use 'set' to remove duplicates, and then sort into order.
(Sidenote: The sort may seem unnecessary, but 1) I believe this is what ensures that $Jup$ corresponds to $J(P)$ listed out in lex-order, which is potentially useful if you're more interested in looking at $J(P)$ rather than just using it as a tool to compute linear extensions and 2) at least for this Python implementation, using 'set' to remove duplicates and then converting back to a list is going to give it to you sorted where you use 'sorted' or 'list', but this makes it more transparent.)
So now we have the list of elements in $Jup$ above $loc[x]$, let's say $K=[a_0,a_1,\ldots,a_r]$. We add new elements to $Jup$, $m+1,m+2,\ldots,m+r+1$ that correspond to the order ideals $a_0,\ldots,a_r$ with $x$ added on. Relations between elements $m$ or less in $Jup$ are unaffected. For relations between the new elements we just added, they will precisely be like the relations in $K$, except reindexed and shifted by $m$. That is to say, if $a_p\lessdot a_q$, then $m+p+1\lessdot m+q+1$.
So now we have to consider relations between elements $m$ or less (corresponding to order ideals with $x$) and elements greater than $m$ corresponding to order ideals that do contain $x$. Luckily, covering relations in $J(P)$ correspond to adding a single element, so the only relations we need to consider are when $x$ is the thing being added. So we add $m+j+1$ to the table of upper covers for $a_j$.
Lastly, we need to update the auxilliary table $loc$. Recall that $loc[y]$ is the principal order ideal of $y$ restricted to elements less than $x$, and now it needs to be the principal order ideal of $y$ restricted to elements less than $x+1$. The only elements affected are elements $y$ that are above $x$ (which we have stored in $up[x]$), as now it's restricted principal order ideal has $x$ added. Since $y\geq x$, before we perform the update, we have that $loc[y]\geq loc[x]$, and so $loc[y]$ was some element of $K$, call it $a_s$. When we update, the new value for $loc[y]$ is the element corresponding to $a_s$ with $x$ added, which is $m+s+1$.
for x in range(n): # We've computed J with P restricted to 0-(x-1), and we're adding x
K = [[loc[x]]]
j=0
while K[j]: #keep adding upper covers of upper covers of upper covers until you can't
K.append([b for a in K[j] for b in Jup[a]])
j+=1
K = sorted(set(item for sublist in K for item in sublist)) #eliminate duplicates, K=[a_0,...,a_r]
for j in range(len(K)): # Adding new element m+1+j corresponding to a_j with x added
i = m+j+1
Jup[i] = [m + K.index(a) + 1 for a in Jup[K[j]]] # any relation between elements of K is still a relation when we
# add x to both elements
Jup[K[j]] = Jup[K[j]] + [i] # relation between a_j and m+j+1 corresponding to adding x
for y in up[x]:
loc[y] = K.index(loc[y]) + m + 1 # update principal order ideal for elements y above x
m += len(K)
So now, all that's left is to compute the number of maximal chains. We described this earlier in terms of building up from the bottom, but since we have a table of up-edges, it will be easier to go from the top-down. As an initial condition, we know there's a unique maximal element (the order ideal consisting of the entire poset). Then we go down through our poset (following the natural labelling in reverse), and we will have that the number of saturated chains from the maximal element to $x$ is the sum of the number of saturated chains from the maximal element to $y$ for all $y$ that cover $x$.
We can actually reuse the same dictionary $Jup$ as we go through, because once we compute the number of saturated chains from the maximal element to $x$, that single number is all the information we need about $x$.
To finish, we realize that $Jup$ also has a unique minimal element (corresponding to the empty order ideal), and the number of saturated chains from the minimal element to the maximal element will be the total number of saturated chains.
Jup[m] = 1
while m>1:
m-=1
ct = 0
for j in Jup[m]:
ct+=Jup[j]
Jup[m]=ct
print(ct)
Now, let's combine this all into a single method, and provide a (hopefully) convincing test example.
def CountLinearExtensions(P):
n = len(P.keys())
up = {key: value[:] for key, value in P.items()}
for i in range(n):
up[n - 1 - i] = sorted(set(up[n - 1 - i] + [item for x in up[n - 1 - i] for item in up[x]]))
Jup = {0: []}
m = 1
loc = [0] * n
for x in range(n): # We've computed J with P restricted to 0-(x-1), and we're adding x
K = [[loc[x]]]
j=0
while K[j]: #keep adding upper covers of upper covers of upper covers until you can't
K.append([b for a in K[j] for b in Jup[a]])
j+=1
K = sorted(set(item for sublist in K for item in sublist)) #eliminate duplicates, K=[a_0,...,a_r]
for j in range(len(K)): # Adding new element m+1+j corresponding to a_j with x added
i = m+j
Jup[i] = [m + K.index(a) for a in Jup[K[j]]] # any relation between elements of K is still a relation when we
# add x to both elements
Jup[K[j]] = Jup[K[j]] + [i] # relation between a_j and m+j corresponding to adding x
for y in up[x]:
loc[y] = K.index(loc[y]) + m # update principal order ideal for elements y above x
m += len(K)
Jup[m-1] = 1
while m>1:
m-=1
ct = 0
for j in Jup[m-1]:
ct+=Jup[j]
Jup[m-1]=ct
return(ct)
The test example included below is the boolean poset corresponding to subsets of 4 elements. To make the structure a little more obvious,write out each number as a 4 digit binary number, and think of it as an indicator vector for a 4 element subset. For example, 6 would be 0110 in binary, so the things that it is covered by would be 1110 and 0111 in binary, or 14 and 7.
P = {0: [8, 1, 2, 4], 1: [9, 3, 5], 2: [10, 3, 6], 3: [11, 7], 4: [12, 5, 6], 5: [13, 7], 6: [14, 7], 7: [15], 8: [9, 10, 12], 9: [11, 13], 10: [11, 14], 11: [15], 12: [13, 14], 13: [15], 14: [15], 15: []}
CountLinearExtensions(P)
Comments About Running Time And Computational Complexity
One thing that's worth mentioning about computing the number of linear extensions of a poset is that we only need to calculate the number of linear extensions on connected components of the Hasse diagram. If you can decompose $P$ into a direct sum of $P_1 \oplus \ldots \oplus P_k$ (ie, there is no relation between $x$ and $y$ unless they are in the same $P_i$, then if $L(P)$ is the number of linear extensions of $P$, and $P$ has $n$ elements, we have \begin{equation} L(P) = \binom{n}{|P_1|,|P_2|,\ldots,|P_k|}\prod_{i=1}^k L(P_k). \end{equation}
As an example, the anti-chain poset (consisting of $n$ elements with no relations) can be thought of as a direct sum of $n$ posets with a single element. Consist with our previous observation that linear extensions in the anti-chain poset correspond to permutations, we get $n!$ without any further work.
Calculating the connected components of the Hasse diagram is a fairly straightfoward graph algorithm, and is computationally trivial compared to the complexity of calculating the number of linear extensions. But in practice, typically one is working with posets that whose Hasse diagrams are already connected. If it decomposed, then those indecomposable pieces would be the things that you would be working with.
Furthermore, there are specific families of posets where the number of linear extensions can be given by more elegant formula. For example, the number of standard Young tableaux of a given partition shape (equivalent to linear extensions of a grid-like poset) can be given via the hook-length formula, and a similar hook-length formula exists for calculating the number of linear extensions of a tree.
More generally, calculating the number of linear extensions of a poset is a long-studied problem, and is known to be a #P-complete problem (Brightwell and Winkler, 1991), so we're not expecting this algorithm to have any sort of provably amazing asymptotic complexity. But in practice, this algorithm performs significantly better than iterating over all linear extensions at the expense of increased memory usage. In SageMath, the algorithm to iterate over all linear extensions is highly optimized and additionally written in Cython, and took about 82s to calculate the number of linear extensions of the above example. By contrast, the cardinality method using this algorithm that I implemented does it in about 1.6ms.