Constructing and sampling directed graphs with given degree sequences

The interactions between the components of complex networks are often directed. Proper modeling of such systems frequently requires the construction of ensembles of digraphs with a given sequence of in- and out-degrees. As the number of simple labeled graphs with a given degree sequence is typically very large even for short sequences, sampling methods are needed for statistical studies. Currently, there are two main classes of methods that generate samples. One of the existing methods first generates a restricted class of graphs, then uses a Markov Chain Monte-Carlo algorithm based on edge swaps to generate other realizations. As the mixing time of this process is still unknown, the independence of the samples is not well controlled. The other class of methods is based on the Configuration Model that may lead to unacceptably many sample rejections due to self-loops and multiple edges. Here we present an algorithm that can directly construct all possible realizations of a given bi-degree sequence by simple digraphs. Our method is rejection free, guarantees the independence of the constructed samples, and provides their weight. The weights can then be used to compute statistical averages of network observables as if they were obtained from uniformly distributed sampling, or from any other chosen distribution.


Introduction and definitions
In network modeling problems [1,2,3,4,5,6,7], one often needs to generate ensembles of graphs obeying a given constraint. A typical constraint is the case when the only information available is the degrees of the nodes, and not the actual connectivity matrix. Note that the node degrees by themselves, that is the degree sequence in general does not determine a graph uniquely: there can be a very large number of graphs having the same degree sequence [8]. Full graph connectivity is uniquely determined by the degree sequence only for a special class of sequences (see Ref. [9] for the case of undirected graphs).
Often, the interest lies in the study of network observables, as determined by the given sequence of degrees, and unbiased by anything else. These can be graph theoretical measures, or properties of processes happening on the network (e.g., spreading processes, such as of opinion or disease). The problem of creating and sampling graphs with a given degree sequence, i.e., degree-based graph construction [10,11], is a well-known and challenging problem that has attracted considerable interest amongst researchers [8,10,11,12,13,14,15,16,17,18,19,20,22,23,24,25,26,27,28,29,30]. There are two main classes of algorithms that are used today to achieve the construction of graphs with given degree sequences. One of them is typically referred to as "switching" or edge-swap based [13,15,16,25,27], while the other one is usually called "matching" or stub-matching based [8,12,14,17,31,32,33,26,34]. Switching methods repeatedly swap the ends of two randomly chosen edges within a Markov Chain Monte-Carlo (MCMC) scheme until a new, quasi-independent, sample is produced. Unfortunately, the mixing time of MCMC schemes for arbitrary sequences is not known in the general case. The other class consists of direct construction methods, which perform pairwise matchings of the half-edges emanating from randomly chosen nodes until all edges are realized. Unfortunately, this method can easily generate multiple edges and self-loops, i.e., edges starting and ending on the same node, after which the sample must be rejected in order to avoid biases [35]. For a comparison of the two classes of methods see Ref. [23].
Recently, a novel degree-based construction [10] and sampling method [11] was introduced for undirected graphs, which has a worst-case scaling of O(N M ), where M is the number of edges (2M is the sum of the degrees, which are given). A similar method was obtained independently in Ref. [34], but that method is less efficient, with a worst-case scaling of O(N 2 M ). Although the algorithm in Ref. [11] is a direct construction method using stub-matchings, it is rejection free, the samples are statistically independent and the algorithm also provides a weight for every realization.
In many systems the interaction between two entities is not mutual but has a direction from one to the other, such as in the cases of human relationships in social networks [36], gene interactions in regulatory networks, trophic interactions in food webs [20,21], etc. Such systems require a representation by directed graphs (digraphs). In fact, undirected graphs can be interpreted as digraphs in which there are two, oppositely directed edges for each connected pair of nodes. Here we present a generalization of the degree-based graph construction problem to directed graphs. Some of the necessary mathematical foundations, laid down in Ref. [30], are here used and expanded to introduce a digraph construction and sampling algorithm. Although the approach follows closely the one introduced by us for the undirected case [11], the generalization is not at all straightforward, and there are significant differences that the directed nature of the links induces. Before we present our algorithm, we introduce some notations, based on Ref [30]. Let us denote by d of non-negative integer pairs, we want to construct a simple directed graph G(V, E) such that node k ∈ V has (d k ) for its in-and out-degrees, respectively, for all k = 1, 2, . . . , N . A simple directed graph is a graph that has no self-loops, nor multiple directed edges in the same direction between two nodes. There can be at most two edges between a pair of nodes, oppositely directed. We call the sequence D a bi-degree sequence (bds for short). When there is a simple digraph with a given bds D for its degrees, we say that the bds is graphical and that the digraph realizes D. Equivalently, we will also talk about "graphicality" as a property. We distinguish realizations as labeled digraphs, and do not deal here with isomorphism questions. That is, if two realizations are identical up to a permutation of their indices, i.e., they are isomorphic, we will still consider them distinctly. In order to avoid isolated nodes, in the following we will assume that d Notice that even if a bds is graphical, not all connection sequences are guaranteed to end up with a simple digraph. For example, Fig. 1d) shows a simple digraph realization of D 5 = {(0, 1), (2, 0), (1, 2), (2, 2)}. However, if we were to place the first four edges as in Fig. 1e), we would break graphicality: from there on, we would not be able to complete the realization of the bds without creating either self-loops or In--stubs Out--stubs Figure 2. The construction of a digraph realizing a given bds proceeds by connecting the out-stubs of the nodes to the in-stubs of other nodes. In this "bipartite" representation the vertical grey bars represents single nodes.
multiple edges. Hence, it is important to find an algorithm that builds digraphs with a given bds. As we will see, this is a challenging problem in itself. An algorithm that builds a digraph from a given bds sequentially connects the out-links of a node to the in-links of others. We can think of these out-and in-links as "out-stubs" and "in-stubs" emanating from a node, that are paired up with the corresponding stubs of other nodes. An intuitive representation of this is shown in Fig. 2. As the graph construction algorithm proceeds, the number of stubs of the nodes decreases. At any time during this process we will call the number of remaining instubs and out-stubs of a node its residual in-and out-degrees, and the corresponding bi the residual bds. Finally, another concept we will need to use in what follows is the notion of normal order [30], which is essentially the lexicographic order on the bds. That is, we say that a bds is in normal order, if for all 1 ≤ j ≤ N −1, we have either d j+1 . Thus, the bds D 6 = {(5, 2), (4,4), (4,3), (2,5), (2,4), (2, 1)} shown in Fig. 2, is arranged in normal order. Once a bds is in normal order, we will use the words 'left' or 'right' to describe the directions towards lower or higher index values in the sequence.
The remainder of this paper is organized as follows: Section 2 introduces the fundamental mathematical notions and algorithmic considerations that are at the basis of our digraph construction algorithm. Section 3 presents the algorithm and its derivation details. Readers interested only in the algorithm itself may skip Subsection 3.1 and proceed to the summary described in the beginning of Section 3 and in Subsection 3.2. Section 4 deals in detail with the digraph sampling problem, provides the derivation of the sample weights and presents a simple example. Section 5 is dedicated to the complexity of the algorithm, and Section 6 concludes the paper.

Mathematical foundations
As seen from the examples above, not all sequences of non-negative integer pairs can be realized by simple digraphs. The sufficient and necessary conditions for the realizability of a bds are given by the "FR" theorem [37,38,39]: and for all 1 ≤ k ≤ N − 1 : Given a bds, we can easily test if it is graphical using this theorem, and thus we will also refer to it as the "FR test". Condition (1) states that both the number of inand out-degrees for all nodes must be no larger than the number of other nodes it could connect to, or receive connections from. Condition (2) is a consequence of the requirement that every out-stub must join an in-stub somewhere else; the sequence D 3 given in one of the above examples is not graphical because it fails this condition. Condition (3) is less intuitive. Its left hand side is the total number of in-stubs that the group of k highest in-degree nodes can receive. Within this group, a node's outstubs can absorb no more of those in-stubs from the same group than its out-degree or k − 1 (it cannot absorb from itself), whichever is smaller (giving the first sum on the rhs of (3)). Outside of this group, a node cannot absorb more of those in-stubs than its out-degree or k, whichever is smaller (the second sum on the rhs of (3)). Hence, the necessity of (3). For the complete proof see Refs. [38,39]. Note that the example sequence D 4 above fails condition 3 for k = 3. The FR test is the directed version of the Erdős-Gallai (EG) theorem (test) for undirected graphs. An important note is that bi-degree sequences are less constraining than undirected ones. The out-stub of a node is always connected to an in-stub of another, not affecting that node's out-stubs, whereas such distinction does not exist for the undirected case. Alternatively, if we disregard for a moment the directionality of the links and consider the degree of the node to be the sum of its in-and out-degrees, then the corresponding graph realizing the bds can have two edges running between the same pair of nodes, whereas this is not allowed in the undirected case.

Algorithmic considerations
The FR theorem only tests for graphicality, but it does not provide an algorithm for constructing the digraph(s) realizing the given bds. At first sight this might not seem an issue. However, the sequence D 5 in Figs. 1d and 1e) reminds us that graphicality can easily be broken by a careless connection of stubs. Clearly, for the purposes of digraph construction, it should not matter which edges we create first, as long as we make sure that every connection made does not break graphicality. In other words, the possibility to create the rest of the edges, so that a simple digraph results in the end, must always be preserved. Thus, the key for the creation of an algorithm that builds simple digraphs realizing a given bds without rejections is in a theorem that allows us to check if we would break graphicality by placing a specific connection. Indeed, such theorems exist, and they will be discussed below. However, interestingly, they require that connections be made from the same node, until all its stubs are used away into edges. That is, assuming that we already made some connections from a given node i, preserving graphicality, these theorems give necessary and sufficient conditions for keeping graphicality by the next connection still involving node i. Simply put, they won't work in general, if we would attempt a new connection from j to k, where j, k = i, while node i still has dangling stubs.
The connections already made from i to some set of nodes X i represent a constraint for the new connections from i, as these novel connections must avoid the set X i . We call such a constraint associated to a node a star constraint on that node. Once all the stubs of node i are connected into edges while preserving graphicality, we obtain a graphical residual sequence D on at most N − 1 nodes. Clearly, the new connections we make from this point on will not be constrained in any way by the connections we made from node i. For the purposes of realizing the sequence D we can just simply remove node i with its fully completed connections, create a realization by a simple graph of D , then, in the end, add back node i with its connections to this graph in order to obtain a realization of D. The comments above hold both for the undirected and directed cases.
One might think of using the EG test for the undirected case and the FR test for the directed case on a residual degree sequence to decide if graphicality was broken after attempting a new connection from the same node. For the undirected case, we have shown in Ref. [11] that the passing of the EG test by the residual sequence is only a necessary condition, if there is already a star constraint on a node. For example, consider the graphical degree sequence d = {6, 5, 5, 3, 3, 2, 1, 1}, and assume that we made connections from node i = 3 to nodes X 3 = {1, 6, 7}. The residual sequence after these connections is d = {5, 5, 2, 3, 3, 1, 0, 1}. It is easy to check that it passes the EG test. However, we will break graphicality with every realization of d , because it will form a double edge with one of the existing connections from node i = 3 to X 3 . Thus, additional considerations have to be made to ensure the graphicality of the residual sequence for the undirected case, as described in [11]. For the directed case here we use the sufficient and necessary conditions for graphicality under star constraints as provided by Theorem 2 below, proven in Ref. [30].
From now on, we will always talk about algorithms that first finish all the outstubs of a node before moving onto another node with non-zero out-degree. In the case of a graphical bds, once all the out-degrees of all the nodes have been connected into directed edges, we are guaranteed to have completed a digraph, because the total number of in-stubs equals the total number of out-stubs, according to property (2).

Theorems on which the algorithm is based
An algorithm that builds graphical realizations of degree sequences of simple undirected graphs is the Havel-Hakimi (HH) algorithm [40,41]: we choose any node with nonzero residual degree, then we connect all its stubs to nodes with the largest residual degrees avoiding self and multiple connections. This process is repeated with other nodes until all stubs of all nodes are used. There is a corresponding version of the HH algorithm for bi-degree sequences as well, introduced first in Ref. [42], then rediscovered independently in Ref. [30], the latter providing an alternative proof. The HH algorithm for bds proceeds as follows: given a normal-ordered bds, choose any node with non-zero residual out-degree, then connect all its out-stubs to nodes with the largest residual in-degrees, without creating multiple edges running in the same direction, nor self-loops. Reorder in normal order the residual sequence and repeat this process until all stubs of all nodes are used. While for any given bds, the HH algorithm will construct a set of digraphs, it cannot construct all possible digraphs realizing the same sequence, as shown in Ref. [30]. For example, the HH algorithm can never result in the digraph shown in Fig. 1c) realizing the example sequence D 2 above. It is easy to see why: there are two kinds of nodes in this example, with bi-degrees (3, 0) and (1,2). The only nodes with non-zero out-degrees are the (1, 2) types. Using the HH algorithm, we would have to connect both out-stubs of such a node to the nodes with the largest in-degrees, that is to the two (3, 0) types. However, the digraph in Fig. 1c) does not have a (1, 2) node being connected to both (3, 0) nodes, yet it realizes the sequence. The limitation of the HH algorithm comes from the fact that it prescribes to connect the out-stub of a node i to an in-stub of the node with the largest residual in-degree that does not yet receive a connection from node i. However, there can be other nodes whose in-stubs can form a connection with an out-stub of i without breaking graphicality. This shows the importance of finding a method able to build not just a realization of a bds, but all the possible realizations of any given bds.
In the remainder, given a residual bds D, we denote by A i D the allowed set of i, i.e., the set of all nodes to which an out-stub of i can be connected without breaking graphicality. Also, let us denote by X i D the set of nodes to which connections were already made from i, thus representing the star constraint at that stage.
The graphicality test under a star constraint on node i is provided as Theorem 2 below. In order to announce it, however, we need to introduce one more definition. Consider a bds D and a given node i with out-degree d where |S| denotes the number of nodes in S, i.e., its size, and for every node j ∈ S, d (i) j > 0. Next, we take D and reduce by unity the in-degrees of all its nodes in S, then reduce by |S| the out-degree of node i. The bds D thus obtained will be called the bds reduced by S about node i from bds D. Equivalently, D is the residual sequence obtained from D after connecting an out-stub from i to an in-stub of every node from S. Theorem 2 (Star-constrained graphicality) Let D be a bds in normal order on N nodes, and let i , be a set of nodes whose in-stubs are forbidden to be connected to the out-stubs of node i (including i). Define L i as the set of the first ("leftmost") d The proof of this theorem is found in Ref. [30]. What this theorem does is to turn a star-constrained graphicality problem for bds D into an unconstrained one on the reduced bds D . The graphicality of D is then easily tested via the FR theorem. The set L i as defined above will be called the leftmost set for node i.
Although announced in its full generality, as X i could be any predefined subset of i , this theorem applies directly to the digraph construction process when X i represents the set of nodes to which connections were already made in previous steps from the same node i, hence forbidding us to make further connections from i to these very same nodes. In this case, the bds D represents the residual sequence D at that stage of the construction process.
As discussed above, in order for us to be able to construct all the simple digraphs that realize a given bds, we need to find the allowed set A i D for the next outstub of i. Clearly, after every connection from the same node i, the residual sequence changes, and along with it the allowed set may change as well. In order to find A i D for the next out-stub of node i, we could just simply attempt connections sequentially to every node with non-zero in-degree not in X i D ∪ {i}, and test for graphicality after each attempt using Th. 2. The set of nodes for which graphicality would have been preserved would form A i D .
However, this would be inefficient and, actually, not needed. In fact, we can exploit a result which states that, if graphicality is broken by a connection, it will be broken by all other connections to the right of the previous one, in the normal order sense. This is expressed in the following: Theorem 3 Let D be a graphical bds in normal order and let X i be a forbidden set for node i, with i ∈ X i . Let j < k be two nodes such that j, k / ∈ X i . If the residual bds D j obtained from D after forming an edge directed from i to j is not graphical, then the bi-degree sequence D k obtained from D by forming a directed edge from i to k is also not graphical.
This theorem follows from the direct contraposition of Lemma 6 in Ref. [30]. Thus, what we need to do is to find efficiently the leftmost node q in the residual sequence in normal order, a connection to which would break graphicality. We will refer to this node q as the leftmost fail-node. All connections to this node and to nodes to its right are guaranteed to break (star-constrained) graphicality, whereas all connections to its left (with the exception of forbidden nodes and self) are guaranteed to preserve the graphical character.
Note that both theorems 2 and 3 are based on the HH theorem for bi-degree sequences. In fact theorem 2 is a generalization of the HH theorem to include star constraints. Also note that, while for the FR theorem only the in-degrees must be ordered non-increasingly, for the HH theorem and hence for both theorems 2 and 3, the bds must be in normal order, as ordering by in-degrees only is not sufficient. This is easily seen from the following example of graphical bds (not in normal order) D 7 = {(2, 0), (2, 1), (0, 1), (0, 2)}. Using the HH theorem, if we do not worry about normal ordering, but just order by in-degree, we could choose to connect the outstub of node (0, 1) to an in-stub of node (2, 0), then the out-stub of node (2, 1) to the remaining in-stub of (2, 0) (connecting to the largest residual allowed residual indegree), after which we have clearly broken graphicality: both out-stubs of (0, 2) now must be connected to the two in-stubs of (2, 1).
We are now ready to present our digraph construction algorithm, which produces random samples from the set of all possible simple digraphs realizing a given bds.

The algorithm
Given a graphical bi-degree sequence D in normal order (initially D = D): 1) Define as work-node the lowest-index node i with non-zero (residual) out-degree. 2) Let X i be the set of forbidden nodes for the work-node, which includes i, nodes with zero in-degrees and nodes to which connections were made from i, previously. In the beginning, X i includes only the work-node and zero in-degree nodes. 3) Find the set of nodes, A i that can be connected to the work-node without breaking graphicality. 4) Choose a node m ∈ A i uniformly at random and connect an out-stub of i to an in-stub of m. 5) After this connection add node m to X i . 6) If node i still has out-stubs, bring the residual sequence in normal order, then repeat the procedure from 3) until all out-stubs of the work node are connected away into edges. 7) If there are other nodes left with out-stubs, reorder the residual degree sequence in normal order, and repeat from 1).
The most involved step of the algorithm is finding the allowed set (step 3)), which is described next.

Finding the allowed set
Let i be the work-node chosen as in 1) and let D denote the normal ordered, residual sequence obtained after having connected some of the out-stubs of i to in-stubs of other nodes, such that graphicality is still preserved. These previous connections from node i form the set of forbidden nodes X i for the next out-stub σ of i. X i also contains the work-node i itself i ∈ X i and all other nodes with zero in-degrees. Let L i be the set of the first (lowest index) d (o) i nodes from D, not in X i . As D is (star constrained) graphical, we can connect σ to any of the nodes in L i without breaking graphicality (due to Theorem 2), hence L i ⊆ A i .
Let m be the last element of L i in the normal ordered bds D and let us "color" (label) red all the non-forbidden nodes, i.e., all the nodes not in X i , to the right of node m. Please note that these color labels are associated with the nodes, defined by their bi-degrees, and not with their indices of location in the sequence. This set of red nodes R i forms the set of candidates for the leftmost fail-node q. All other nodes are colored (labelled) black. To find the leftmost fail-node we could simply connect out-stub σ to an in-stub of a red node , add the new connection temporarily to the set of forbidden nodes, bring the new residual sequence into normal order, then test for graphicality using Theorem 2. This procedure could then be repeated sequentially, with going over all the red nodes from left to right, until graphicality would fail for the first time at = q. However, the considerations in the following paragraphs allow us define a better method.
For the sake of argument let us perform the sequential testing as explained above. It would imply the following steps for a given red node : (a) Reduce the out-degree at the work-node i and the in-degree at by unity, that is d resulting in a new residual bds D .
(b) Bring D into normal order (required by Theorem 2). Note that is the only node whose in-degree has changed and only the work-node had its out-degree changed (its in-degree was not affected). Thus, when bringing D into normal order, the relative positioning of all the other nodes does not change. The work-node might have shifted to the right to a new position i within the block of nodes with the same in-degree (i ≥ i), and the red node's new position might have also moved to the right in the normal ordered sequence ( ≥ ). (c) Add to the forbidden set for the work-node. (d) Now, as required by Theorem 2, reduce by unity the in-degrees of the nodes in the left-most adjacency set L i , and reduce the out-degree of the work-node i to zero. This results in the new sequence D . (e) Order the bds D by in-degrees, non-increasingly. (f) Apply the FR theorem to test for graphicality.
Thus, whether the connection of the work-node i to breaks graphicality, ultimately depends on whether the residual bds D fails (or passes) the FR test. However, as we noted before, for the FR test we do not need to have the bds D in normal order, we only need to have it ordered non-increasingly by the in-degrees. Additionally, observe that in step (d) the reduction of the in-degrees always happens on the same set of nodes, independently of the red node , that is the left-most adjacency set L i is the same for all . Thus, in this particular case of Theorem 2's application, ultimately we do not need to bring D into normal order (step (b)), only non-increasingly by in-degrees, which would be done anyway in step (e). That means we can just skip step (b), we do not need to move around any of the nodes at that stage. Thus, the only difference between the sequences D for different -s is at the position of this node after the reordering in (e), with respect to the rest of the sequence.
These observations suggest that we should define a bds D obtained from the bds D by reducing by unity the in-degrees of all nodes in the set L i \ {m} and by d i − 1 the out-degree of the work-node i, leaving only one out-stub (out-stub σ) at i. Clearly, the bds D is graphical (connecting out-stub σ to an in-stub of node m surely preserves graphicality, by Theorem 2). Let us now order D non-increasingly by its in-degrees, in a specific way, described as follows. Shift only the reduced in-degree nodes in D to the right with respect to the rest of the sequence such as to restore non-increasing ordering by the in-degrees (if needed). Since only the in-degrees of the nodes in the set L i \ {m} have been reduced, keep the relative ordering of all other nodes in D exactly the same as in D. Thus the relative ordering of the red nodes and of the work node have been preserved as well. Let us denote the new location of the work node in D by j (j ≤ i). Connecting now σ to an in-stub of a red node in this sequence will produce the same set of residual bi-degrees as in step (d) above. To be able apply the FR theorem, then all we need to do is to shift to the right node in the sequence (if needed) to make sure that it is non-increasingly ordered by in-degrees. Since only the in-degree at was modified (reduced by unity), this reordering is very simple: if x denotes the location of the last node of the block of nodes with the same in-degree as node in D (x ≥ ), then we simply swap the node at with the node at x after the reduction of the in-degree at . Let us denote the obtained sequence by D . Clearly, it is non-increasingly ordered by in-degrees, and thus we can apply the FR theorem to see if it is graphical. Note: it could happen that x = j (e.g., there are many nodes with zero out-degree but larger in-degree than the work-node as defined in 1)), however, the steps below can be applied just the same.
Next, we show how to identify the leftmost red fail-node q by investigating how the inequalities in (3) break down. Since D is graphical, we have for all 1 ≤ k ≤ n − 1 (n is the last element of D ) that L (k) ≤ R (k), where L and R are the left hand side (lhs) and the right hand side (rhs) of inequalities (3) written for D : Let us denote by L (k) and R (k) the lhs and rhs of the inequality (3) corresponding to D . Since the rhs of (3) involves only out-degrees, and we only reduced the outdegree of the work-node from 1 to 0, we will always have R (k) = R (k) − 1, except when k = 1 and the work-node is at j = 1, in which case R (1) = R (1). However, in this case, L (1) = L (1), because only the in-degree of j = 1 appears, which does not get changed. Thus, since L (1) ≤ R (1) (D is graphical), graphicality cannot be broken at k = 1 when j = 1. Let us now consider that the work-node is still at position j = 1, but k > 1. For 1 < k < x, the in-degrees in D are the same as those in D , hence L (k) = L (k). For k ≥ x, however, we have L (k) = L (k) − 1. Now consider j > 1. For 1 ≤ k < x, we have L (k) = L (k) and for k ≥ x, L (k) = L (k) − 1. The following summarizes the relationships above: Since L (k) ≤ R (k) for all k, graphicality for D can only be broken (that is to have L > R for some k), if L (k) = R (k), namely in cases (A.2) and (B.1) above.
Observe that L (k) and R (k) are computed from D , hence they are independent from or x. This gives us the following simple procedure for finding the leftmost fail-node, if it exists. Starting from k = 2 for j = 1, and k = 1 for j > 1, find the smallest k 0 for which L (k 0 ) = R (k 0 ). If no such k 0 exists, then there are no fail-nodes and all non-forbidden nodes are to be included in the allowed set. If there is such a k 0 , the first red node q to the right of k 0 (q ≥ k 0 + 1) is the leftmost fail-node of D , which when identified in the original bds D will give the leftmost fail-node q. All non-forbidden nodes to the left of q are to be included in the allowed set.

Summary for finding the allowed set
What we discussed in detail in the previous subsection corresponds to step (3) of the main algorithm described in the beginning of Section 3. Given the normal-ordered bds D at the end of step 2) of the main algorithm: nodes not in X i . (3.2) Identify the "red" set R i as those nodes that are neither in L i nor in X i . Note, the color label is associated with the node, not its index. (3.3) Build D as follows: where m is the last node in L i . (3.4) Shift nodes from L i \ {m} to the right in the sequence (and only these) such as to restore ordering non-increasingly by in-degrees (if needed), preserving the color labels of the nodes in the process. The work-node may have shifted to a new location j after this step. This is the updated sequence D .
(3.5) Starting from k = 1 if j = 1 or from k = 2 if j = 1, find k 0 as the smallest k such that L (k) = R (k), where L (k) and R (k) are computed from the reordered (after step (3.4)) D using (4) and (5). If there is no such k 0 , then the allowed set A i is all the nodes in D except nodes from the forbidden set X i . (3.6) Otherwise find the leftmost red node q in the updated bds D to the right of k 0 , that is with q > k 0 . Then the corresponding node q in D , will be the leftmost fail node. Note that q is the new position of the node at q in D after the reordering in (3.4). (3.7) The allowed set A i is formed by all nodes in D not in X i , and to the left of q.

The sampling problem
The algorithm generates an independent sample digraph every time it runs, without restarts or rejections, and it guarantees that every possible realization of a graphical bds by simple digraphs can be generated with a non-zero probability. However, the algorithm realizes the digraphs with non-uniform probability. Nevertheless, knowing the relative probability for every digraph's occurrence allows us to calculate network observable averages as if they were obtained from a uniform sampling. In particular, the following expression, which is a well-known result in biased sampling [43,44], provides these averages as: where Q is an observable measured from the samples s j generated by an algorithm. The w(s j ) sample weight is the inverse of the relative probability of the occurrence of s j and M is the number of the samples generated. In Subsection 4.1 we give a detailed derivation of this formula, specialized to our graph construction problem. The weights of the samples generated by our algorithm are given by where i runs over all the nodes with non-zero out-degree as they are picked by the algorithm to become work-nodes, and k i (j) = |A i (j)| is the size of the allowed sets A i (j) just before connecting the j-th out-stub of i. Note that w ≥ 1 since there always exists at least one digraph realizing the bds. Subsection 4.2 gives a derivation of (7).

Biased sampling over classes
Our algorithm sequentially connects all stubs from a series of work nodes and finishes with a simple, labeled digraph. This process can be uniquely described by a path of connection sequences. Having chosen a work node i 1 for the first time, it determines the allowed set A i1 . We next choose uniformly at random a node j 1 (i 1 ) ∈ A i1 and connect a stub of i 1 to a stub at j 1 (i 1 ). We could have chosen j 1 (i 1 ) following any other criterion, but in that case the expression (7) of the weights would have to be modified accordingly. After this connection we recompute the new allowed set A j1 (i 1 ), then connect another stub of i 1 , and so on until all the stubs have been used up at i 1 . Let us denote by s such a path of connection sequences: whered (o) i denotes the residual out-degree of node i. A path s uniquely defines the digraph G(s) created, as the collection of all connections in (8) forms the edge set of the created graph G(s). However, several paths may lead to the same digraph. Also note that the order of the connections in (8) matters in the calculation of the weight, as the corresponding allowed sets in general depend on history of connections up to that point. For a finite bi-degree sequence the number of distinct samples (paths) is also finite. Let us denote this set of paths by: Let us now assume that we built with our algorithm a sequence of samples s 1 , s 2 , . . . , s M , and that the sample number M is large enough for us to see all elements of Π sufficiently many times. Given some path s we compute a quantity Q(s), and we are interested in calculating the average of Q over path ensembles. In our case Q is defined on the final graph itself Q(s) = Q(G(s)), but for now we will not consider that, explicitly. If we were just simply computing the average of Q over the set of samples, we would obtain an average biased by the way the algorithm builds the paths from Π: where M k is the number of times we have seen path π k appear in the sequence of samples. Clearly, is the probability by which path π k is generated via the algorithm. We now assume that we can compute analytically the path probabilities ρ k , from knowing how the algorithm works. Instead of (9) we want to compute the average as if it was measured over the uniform ensemble of paths, that is: If we form: , we have lim M →∞ Q bp = Q up , due to (10). Thus, the weighted average (12) should be used in order to obtain averages according to uniform sampling in the M 1 limit. Let us assume that there is an equivalence relation "∼" between paths, hence inducing a partitioning of Π into K equivalence classes: Π = C 1 ∪ . . . ∪ C K , where C = π k 1 , . . . , π k µ . The size of class C is denoted by µ = |C |. We have K =1 µ = P . Alternatively, for some given path π, we will denote by C(π) the equivalence class of π and by µ(π) = |C(π)| its size. Let us also assume that if s, r ∈ C , that is s ∼ r, then Q(s) = Q(r). For example, in our case distinct paths may lead to the same digraph. We introduce the equivalence relation "∼" and say that two paths s and r are equivalent, s ∼ r if they produce the same labeled digraph, G(s) = G(r). Clearly, if Q depends only on the constructed graph, i.e., Q(π) = Q (G(π)) for all π ∈ Π, then Q(s) = Q(r) whenever s ∼ r.
Our goal is to obtain the average of Q uniformly over the equivalence classes, that is: where we chose to write the first element of C in the argument of Q, but of course, any other element could have been chosen from the same class, as Q is constant within a class. In general, (12) will not produce Q uc , but a sum weighted by class sizes. Instead, let us consider: It is then easy to see that: In order for (14) to be useful in practice, one has to be able to compute the size of the equivalence class µ(s) from seeing s and knowing how the algorithm works. Fortunately this is possible in our case, as shown next.

Computing the weights
First, let us note that when connecting the out-stubs of a work-node we are not affecting the out-stubs of any other nodes, but only in-stubs. Hence, all nodes with non-zero out-degrees will eventually be picked as work-nodes by the algorithm. Since normal ordering is first by in-degrees, the order in which nodes will become worknodes depends on the sequence of connections. Let us now calculate the probability of the path s in (8). Given a residual sequence, the work-node i 1 is uniquely determined by the algorithm as described before. Since the next connection is picked uniformly at random, the probability of the link from i 1 to j 1 (i 1 ) ∈ A i1 (j 1 ) is |A i1 (j 1 )| −1 . Let k i (j) = |A i (j)| denote the number of nodes in A i (j). Then, it is easy to see that the probability of a path s is given by: where i 1 , i 2 , . . . , denote the work-nodes in the order in which they are picked by the algorithm. This expression can be computed readily in a computer as the algorithm progresses. In order for us to use (14) it seems that we would need also to obtain the size µ(s) of the class to which path s belongs. Clearly, two different paths s and s will result in the same graph (s ∼ s ) if and only if the sequence of connections in one path is a permutation of the connections in the other path. Hence, the class size µ(s) is nothing but the number of permutations of the connections, which is the same for all paths, that is, all classes have the same size µ. Since all connections are made from a node first before moving on to another, we have µ = we actually don't need to use this number: one can simply multiply by µ both the numerator and the denominator of (14) to obtain (6-7).

A simple example
In this subsection we illustrate the algorithm on a simple sequence: Let us now consider the Pearson coefficient r of degree-degree correlations, or the assortativity coefficient defined for directed graphs [45] as our network observable Q = r. For each one of the 11 graphical realizations of D 8 , r can be calculated exactly, as can the uniform average over this ensemble, obtaining r theo = −0.040506. We will refer to r theo as the "theoretical value". We then let our algorithm run on this sequence to produce M samples and using (6)(7) to obtain the corresponding coefficient r M . Fig 3a) shows a few runs with different seeds and their convergence to the theoretical value. Fig 3b) shows the standard deviation ([ r M − r theo ] 2 ) 1/2 where the overline denotes an ensemble average over runs.

Complexity of the algorithm
To determine the theoretical upper bound for the complexity of the algorithm, that is the worst-case complexity, notice that there are only three steps in the algorithm that require more than O (1) computational operations, or steps, to complete. First, after each connection is placed, one must bring the residual sequence into normal order, steps 6) or 7). To accomplish this, both the work-node i and the target node m will have to move to the right, but the relative positions of all other nodes will remain unchanged. In other words, if we were to remove nodes i and m, the rest of the bds would already be sorted. Thus, in order to complete these steps, one only has to find the new positions of nodes i and m and insert them into the already sorted bds. Therefore, the complexity of either one of step 6) and step 7) is simply where N is the number of the nodes in the sequence being ordered.
Second, the allowed set A must be built before placing each connection (step 3). Following the summary of this step, given in Subsection 3.2, notice that steps 3.1 to 3.4 can be all finished during a single scan of the residual bds. This is clearly so for the creation of the leftmost set L i and for setting the "red" color labels (or flags) (steps (3.1) and (3.4)). Concerning the ordering of the bds D , it is possible to create it already sorted by simply scanning the bds D while keeping track of the in-degree d of the nodes currently being copied and the index a in D of the first node with that in-degree. Then, because D is in normal order, the only possibility for a node in D to break the order is if its in-degree equals d + 1. In this case, it can be simply swapped with the node at a, because, as argued in Subsection 3.1, the mechanism to build the allowed set is entirely based on the FR theorem, which does not require the bds to be in normal order, but to be simply ordered non-increasingly by its in-degrees. Thus, steps (3.1) to (3.4) can be completed in O (N ) steps.
Third, the computation of the sums L , R and their comparison must be conducted, which is the same step as (3) in an FR test. To determine the complexity of an FR test note that computing the repeated sums for each one of the inequalities (3) is quite inefficient. Instead, below we derive recurrence relations that allow us to complete the FR test in a linear, O (N ) number of steps.
The steps of the main algorithm are done sequentially, and thus can all be completed in a total of O (N ) steps. They must, however, be repeated for each edge in the digraph. Thus, the maximum complexity of the algorithm is i is the number of edges. Since O (M ) ≤ O N 2 the maximum complexity of the algorithm is O(N 3 ). It is important to note though, that for a given bds the complexity of the algorithm can be substantially smaller, similar to the case for our undirected graph sampling algorithm [11].

The Fulkerson-Ryser test revisited
The most complex part of the Fulkerson-Ryser test is to compute the lhs and the rhs of inequalities (3), which we rewrite here for the sake of readability: Our goal is to find recursion relations for L(k) and R(k). For the lhs the relation is simply For the rhs, first note that one can write it as where g (o) i (k) is the family of integer sequences defined as , that is, the number of indices i for which g i (k) = p. Then, from (16) follows that hence where we used the fact that N p=0 G k (p) = N . Furthermore, let us introduce the following notations: Then, after some simple manipulations, from (17) it follows that Finally, notice that ∆G k (p) = δ p, Substituting it into (19), we obtain: Thus, we have turned the problem of finding a recursion relation for R(k) into the problem of findingG k (k). To solve this, first note that withG 1 (1) = G 1 (0) + G 1 (1). The above equation constitutes a recursion relation for G k (q). Such a relation can be rewritten as Observe that S (k) and G 1 (k) can be easily computed while scanning the bds, and then calculating L(k) and R(k) for each k requires a single operation. Thus, the entire FR test can be completed in O (N ) steps.  Figure 4. Probability distribution p of the logarithm of weights for an ensemble of bi-degree sequences on N = 100 nodes. The in-degrees were drawn from a normalized power-law distribution ∼ d −γ in with γ = 3 and the out-degrees were drawn from a Poisson distribution e −λ λ d out /dout!, with the same average as the average in-degree, λ = d in . The black circles are the simulation data and the red continuous line is a Gaussian fit.

Discussion
In summary, we have developed a graph construction and sampling algorithm to construct simple directed graphs realizing a given sequence of in-and out-degrees. Such constructions are needed in practical modeling situations, ranging from epidemics and sociology through food-webs to transcriptional regulatory networks, where we are interested in learning about the statistical properties of the network observables as determined only by the bi-degree sequence and nothing else.
Unlike existing algorithms such as the Configuration Model, which is affected by uncontrolled biases and unacceptably long running times except for a very restricted class of sequences, our algorithm is rejection-free. Also, it guarantees the independence of the produced samples, unlike MCMC methods, which have unknown mixing times. While its mathematical underpinnings are nontrivial, the algorithm itself is straightforward to implement. In principle, our approach can be extended to include more complex constraints, such as a given sequence of motifs frequencies, but we have only concentrated on degree sequences since they are, arguably, the most fundamental of constraints. The algorithm can also be used to sample from given in-and outdegree distributions, not just sequences: given such distributions, one first samples a graphical bds from these, then one applies our algorithm to generate digraphs. In this case, however, the sample weights (7) must be modified to reflect the probability of the occurrence of the given graphical bds when drawn from the distributions.
Just as in the case of undirected graphs, we can expect the distributions of the weights for large graphs to be log-normal, as shown in Ref. [11]. As an example, figure 4 shows the distribution for bi-degree sequences in which the in-degrees follow a power law with exponent γ = 3 and the out-degrees a Poisson distribution whose mean matches the average in-degree. Indeed, the distribution of the weight logarithms is well approximated by a Gaussian. Similarly the undirected case, we find for all the examples we studied numerically, that the standard deviation σ of the distributions of weight logarithms grows slower than the mean m with the number of nodes N ; see figure 5 showing the scaling of m and σ for bi-degree sequences in which both indegrees and out-degrees follow a power law distribution with exponent γ = 3. Thus, we may expect that typically, in the N → ∞ limit, the rescaled weight distribution becomes a delta function, making the sampling asymptotically uniform. Bounds on the complexity of the algorithm could easily be obtained by inspecting the algorithm, showing a maximum complexity on the order of O(N M ) where M is the total number of edges, M = N i=1d (o) i . In developing our results, we also provided an efficient way of implementing the Fulkerson-Ryser test, whose scope of application goes beyond our present algorithm, as it can be used in any context to determine whether a bi-degree sequence is graphical. copyright notation here on.