{'fraction_non_alphanumeric': 0.034858533008964576, 'fraction_numerical': 0.012463758456360182, 'mean_word_length': 4.666163141993958, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Derivative structure enumeration using binary decision diagram

A derivative structure is a nonequivalent substitutional atomic configuration derived from a given primitive cell. The enumeration of derivative structures plays an essential role in searching for the ground states in multicomponent systems. However, it is computationally hard to enumerate derivative structures if the number of derivative structures of a target system becomes huge. In the present study, we introduce a novel compact data structure of the zero-suppressed binary decision diagram (ZDD) to enumerate derivative structures much more efficiently. The present study shows its simple applications to the enumeration of structures derived from the face-centered cubic and hexagonal close-packed lattices in binary, ternary, and quaternary systems. The present ZDD-based procedure should significantly contribute not only to various computational approaches based on the derivative structures but also to a wide range of combinatorial issues in physics and materials science.


I. INTRODUCTION
Structure enumeration has played an essential role in performing crystal structure prediction and in understanding crystal structures. In general, structure enumeration requires a given policy that restricts an entire set of structures in a continuous configuration space to a discrete set of structures. Such a policy is the atomic substitution of a given structure, and nonequivalent substitutional structures are called "derivative structures" [1]. Although the convex hull in a configurational space [2,3] and an efficient algorithm specially developed to find the ground states in multicomponent systems [4] have recently been found to be applicable when the cluster expansion method [2,5,6] can describe the configurational energy accurately, a set of derivative structures has been commonly used to search for the ground states in multicomponent systems (e.g., Refs. [7][8][9]. A set of derivative structures itself is also of interest from the viewpoint of crystal chemistry because many existing crystal structures of not only intermetallic alloys but also ionic compounds have been interpreted as derivative structures [10,11]. The well-known Pólya counting theorem [12] has a long history of being used to count the number of nonequivalent molecule structures [13] and the number of derivative structures [14], because they can be regarded as graph coloring problems of assigning colors to graph vertices under a given set of permutations. Recently, Hart and Forcade proposed an efficient procedure to enumerate the derivative structures themselves, not only their total number [15,16]. Their procedure is based on the enumeration of nonequivalent lattices represented by the Hermite normal form (HNF) and the enumeration of labelings using a finitely generated Abelian group given by the Smith normal form (SNF). Moreover, faster algorithms have also been reported recently [17,18].
These procedures are practically sufficient to enumerate derivative structures required to determine the ground-state structures in a binary alloy with a simple lattice such as a face-centered cubic (fcc) one, because each of many intermetallic compounds has a primitive cell composed of up to 24 atoms [19]. On the other hand, the possible size of periodicity in derivative structures is very restrictive in a binary system with a small number of symmetry operations and a system with three or more components. In this study, we propose a much more efficient method of enumerating derivative structures. We employ a compact data structure developed in the algorithm theory for representing a set of combinations. In particular, we use the zero-suppressed binary decision diagram (ZDD) [20], which has been used to enumerate constraint subgraphs of a given graph, such as all paths between two given vertices (s-t paths) [21,22] and spanning trees [23]. For example, ZDD succeeded in enumerating s-t paths for the 27 × 27 grid graph, the total number of which reaches as many as approximately 10 163 [24,25]. The use of ZDD should enable us to significantly increase the possible size of periodicity in derivative structures.
This paper is organized as follows. Section II introduces the terminology for representing a derivative structure mathematically. Section III shows fundamental ideas to eliminate equivalent structures among all possible substitutional ones, following the works of Hart and Forcade [15,16]. Section IV introduces ZDD and demonstrates how to apply ZDD to the derivative structure enumeration. The present method employs an isomorphism-eliminated ZDD proposed in Ref. 26, which can be applied to the enumeration of binary derivative structures in a straightforward manner. In addition, we propose a generalization of the isomorphism-eliminated ZDD to extend the scope of application including the derivative structure enumeration in multicomponent systems. Finally, Section V shows the application of the present ZDD-based method to the enumeration of binary, ternary, and quaternary derivative structures from the fcc and hexagonal close-packed (hcp) primitive cells.

II. TERMINOLOGY TO REPRESENT DERIVATIVE STRUCTURES
In this section, we define the terminology required to represent a derivative structure mathematically. We mainly follow the convention of crystallography to define them [27]. Figure 1 illustrates the terminology in a twodimensional example, and their mathematical definitions will be given for three dimensions in the following.

A. Parent lattice and primitive cell
Given a set of basis vectors A = (a 1 , a 2 , a 3 ), the parent lattice L is defined as where Z denotes the set of integers. A primitive cell is also defined as the pair of the set of basis vectors A and the set of sites whose positions are described by point coordinates D, as shown in Fig. 1 (a). We adopt the convention where A is used as the basis vectors to describe the point coordinates. Therefore, the values of the point coordinates in the primitive cell range from zero to one. All of the point coordinates in its crystallographic pattern are expressed as {d + | d ∈ D, ∈ Z 3 }.

B. Sublattice and supercell
A sublattice is a subset of the parent lattice L obtained by removing some lattice points from the parent lattice L [27] [28]. A set of basis vectors of the sublattice is identified with the transformation matrix M such that the original set of basis vectors A is transformed into a new set of basis vectors AM. Therefore, the sublattice L M is the set of lattice points expressed as We refer to the determinant of M, det M, as the index of the sublattice L M . The index is identical to the number The shaded region represents a supercell. Since the index of sublattice (b) is two, the number of sites in the supercell is twice that in the primitive cell. We denote the point coordinates of the sites in the supercell as D M = {d1, d2, d3, d4}. We use matrix A as the basis vectors to describe the point coordinates of the sites in the supercell D M . (d) A binary derivative structure with supercell (c). In our label notation, integers 0 and 1 represent the blue and yellow atoms, respectively. Therefore, labeling c = (0, 1, 1, 0) indicates that the blue atoms occupy sites 1 and 4 and the yellow atoms occupy sites 2 and 3.

C. Derivative structure and labeling
A k-ary derivative structure is defined as a nonequivalent structure in which one of k atomic species occupies every site of a supercell identified with AM and D M . Therefore, a derivative structure can be equivalently regarded as the labeling of the sites for the supercell. The labeling can be expressed using k integers where c i denotes the label of site i. Each label indicates one of the atomic species. Here we follow the labeling in descending order used in Ref. 26. Figure 1 (d) illustrates a derivative structure in the two-dimensional example.

III. EQUIVALENT STRUCTURE ELIMINATION
Since the present ZDD-based method follows two fundamental ideas to eliminate equivalent structures, which were used in the successive works of Hart and Forcade [15,16], we summarize them in this section. One idea is the equivalent sublattice elimination using HNF (Section III A). The other is the equivalent labeling elimination using the structure of the finitely generated Abelian group of a given sublattice (Section III B). In Section III C, we demonstrate a two-dimensional example of the equivalent labeling elimination for a given sublattice.

A. Equivalent sublattice elimination
As described in the previous section, a sublattice is obtained by transforming basis vectors using the transformation matrix M. An infinite number of integer transformation matrices are possible even for a given index, which is closely related to the arbitrariness for choosing basis vectors of a given lattice. Fortunately, however, the number of nonequivalent sublattices is finite for a given index. We can enumerate a complete set of nonequivalent sublattices for the index.
Let U be a three-dimensional square unimodular matrix, where all elements are integers and det U = ±1. It is well known that matrices M and MU are equivalent in terms of lattice transformation [30]. This means that they derive the same sublattice expressed as although they give different sets of basis vectors spanning the sublattice. Their representative can be the canonical form called the Hermite normal form (HNF) [30]. Any transformation matrix M can be converted to a unique form of the lower-triangular integer matrix, HNF, by multiplying the unimodular matrix U from the right satisfying the relationship where a > 0, 0 ≤ b < c, 0 ≤ d < f , and 0 ≤ e < f . The requirement that diagonal elements a, c, and f are all positive eliminates equivalent basis vectors obtained by inversion. Also, the addition of a basis vector to another one or the subtraction of a basis vector from another one does not change the lattice itself. Thus, we can choose remainders of f as d and e, and a remainder of c as b.
To enumerate nonequivalent sublattices for a given index, therefore, it is sufficient only to enumerate HNFs whose determinant is the index. HNFs are easily enumerated by brute force. The product of diagonal elements a, c, and f should be equal to the index and the diagonal elements should be divisors of the index. For each set of diagonal elements {a, c, f }, all combinations of non-diagonal elements b, d, and e satisfying the inequalities can be generated. In what follows, we consider only transformation matrices in the lower-triangle HNF.
Then, we eliminate equivalent sublattices among the enumerated ones according to the symmetry of the crystallographic pattern of the primitive cell, as performed in the works of Hart and Forcade [15,16]. We denote a symmetry operation in the space group of the primitive cell [31]

B. Equivalent labeling elimination
Given the sublattice L M , we enumerate nonequivalent labelings of the sites in the supercell with the sublattice L M . To eliminate equivalent labelings from all possible labelings, we introduce a permutation representation for symmetry operations of the supercell. Considering symmetry operations in the space group of the primitive cell that leaves the supercell unchanged, they form a subgroup of the space group of the primitive cell H M . If the symmetry operation g ∈ H M moves site i to site j in the supercell, we describe operation g as a permutation, σ g (i) = j. By applying this rule to all symmetry operations in H M , we obtain the permutation group Σ M mapped from the space group H M as which is homomorphic to the space group H M . We then define a permutation of labeling c. Given permutation σ ∈ Σ M , permuted labeling σ(c) is expressed as We refer to labelings c and σ(c) as equivalent labelings for the supercell with the permutation group Σ M . Therefore, the set of labelings equivalent to c is given by the orbit of c as A representative of the labeling orbit can be defined as the maximum labeling in the lexicographical order, as described in Ref. 17. For example, labeling c 1 = (0, 1, 1, 0) is larger in the lexicographical order than labeling c 2 = (0, 0, 1, 1). We denote that c 1 is larger than or equal to c 2 by c 1 c 2 . Finally, nonequivalent labelings for the supercell with the sublattice L M can be written as (9) Moreover, Pólya's counting theorem is applicable only for determining the size of C M,k or counting the number of nonequivalent labelings [12,13] (see Appendix A).
The set of nonequivalent labelings depends only on the permutation group. This means that C M,k and C M ,k for M = M having a bijection or one-to-one correspondence are isomorphic if they have isomorphic permutation groups, i.e., Σ M ∼ = Σ M . Therefore, it is sufficient to enumerate nonequivalent labelings only for a complete set of non-isomorphic permutation groups.
Note that the set of nonequivalent labelings given by Eqn. (9) contains superperiodic labelings that can be expressed by a smaller supercell and labelings with less than k atomic species denoted as incomplete labelings in Ref. 15.

C. Two-dimensional example
We demonstrate an example of enumerating derivative structures from a two-dimensional primitive cell. The primitive cell is composed of a square lattice and a site at the origin. The basis vectors of the primitive cell and the point coordinates of the site are expressed as The space group type of its crystallographic pattern is p4mm [27]. We hereafter consider the enumeration of derivative structures for the transformation HNF matrix of where the first and second rows represent the original sequence and the permuted sequence, respectively. As described in Appendix A, the number of nonequivalent labelings with k atomic species is easily obtained from Pólya's counting theorem as In the binary case (k = 2), Pólya's counting theorem indicates that there are six nonequivalent labelings for the present supercell. Although the set of all possible labelings is given by {c | c ∈ {0, 1} 4 }, the six nonequivalent labelings are the maximum labelings in the lexicographical order, expressed as They are identical to the nonequivalent labelings shown in Figs. 2 (b)-(g).

IV. DECISION DIAGRAM
In this section, we demonstrate a ZDD-based method to enumerate nonequivalent labelings much more efficiently. In Section IV A, we describe the binary decision tree representing a family of subsets from a finite number of elements. Then we introduce a compact form of the binary decision tree, ZDD, in Section IV B. Finally, we propose a procedure to construct a ZDD representing a set of binary nonequivalent labelings (Section IV C 1) and a procedure to construct a ZDD representing a set of multicomponent nonequivalent labelings (Section IV C 2).

A. Binary decision tree
A binary decision tree [33][34][35] represents a family of subsets generated from n elements satisfying given conditions. For example, subsets containing exactly two elements from among {a, b, c} are {a, b}, {a, c}, and {b, c}. By fixing the order of the elements, we express the family of the subsets as a binary decision tree. The first panel of Fig. 3 shows the binary decision tree representing the family of the subsets, S = {{a, b}, {a, c}, {b, c}}. The binary decision tree comprises terminal nodes, nonterminal nodes, and directed edges. Each non-terminal  = (c4, c3, c2, c1). The blue contour atoms and yellow atoms correspond to labels 0 and 1, respectively. (h) Binary decision tree representing nonequivalent labelings (b)-(g). The non-terminal node ci corresponds to label ci at site i in the supercell. The broken and solid arrows from the non-terminal node ci indicate 0 and 1 assigned to label ci, respectively. The square nodes labeled by 1 and 0 indicate that the labeling is the maximum or not the maximum, respectively, among its equivalent labelings in the lexicographic order. (i) Irreducible ZDD derived from the binary decision tree (h).
node has two kinds of outgoing edges, the 1-edge and the 0-edge. They respectively indicate whether or not a subset includes the element corresponding to the nonterminal node. Therefore, a path from the root node to a terminal node represents a subset. Then, the binary value of the terminal node called the 1-terminal or 0terminal node respectively indicates whether or not the family of subsets S contains the corresponding subset. Therefore, each of three paths reaching the 1-terminal nodes, called 1-paths, corresponds to each subset in S = {{a, b}, {a, c}, {b, c}}.
A binary decision tree also similarly represents a set of nonequivalent labelings. Figure 2 (h) shows the binary decision tree of the set of nonequivalent labelings in the two-dimensional example shown in Figs. 2 (b)-(g). The non-terminal node c i has the 1-edge and the 0-edge, indicating that the corresponding label is assigned as c i = 1 and c i = 0, respectively. The six 1-paths are identical to the nonequivalent labelings.

B. Zero-suppressed binary decision diagram
A binary decision diagram (BDD) is a canonical representation for a Boolean function [36,37], derived by reducing a binary decision tree to a directed acyclic graph. A ZDD is a variant of the BDD [20,38] and specially designed for representing sets of combinations. ZDDs are more efficient than BDDs for representing a family of sparse subsets [20]. A ZDD is derived by reducing a binary decision tree on the basis of the following two reduction rules, as schematically illustrated in Fig. 4. (1) All nodes whose 1-edge directly points to the 0-terminal nodes are deleted, and (2) all equivalent nodes having the same child nodes and the same variable are shared. eliminated, as can be seen in the second panel of Fig. 3. Then, redundant terminal nodes are combined into single 1-terminal and 0-terminal nodes, as shown in the third panel of Fig. 3. Finally, two non-terminal nodes c with the same child terminal nodes are shared, as shown in the fourth panel of Fig. 3. The obtained irreducible ZDD is more compact than the binary decision tree in terms of the number of nodes. Moreover, the irreducible ZDD is canonical and independent of the order of the reduction processes for a given order of variables.
Each 1-path in a ZDD corresponds to a solution. The following three 1-paths are identical to the family of subsets S = {{a, b}, {a, c}, {b, c}}.
Note that all paths from the root node to the terminal nodes are represented with a common order of variables. Figure 2 (i) shows the irreducible ZDD for the set of nonequivalent labelings constructed from the binary de-cision tree shown in Fig. 2 (h). The following six 1-paths in the ZDD represent the six nonequivalent labelings. They are identical to the set of nonequivalent labelings given by Eqn. (14). In general, the number of paths is calculated by dynamic programming in a computational time proportional to the number of ZDD nodes. This is typically a fast way to count the number of paths because the number of nodes of a ZDD is much smaller than that of paths in most cases.

Binary system
Here, we reformulate the definition of the set of nonequivalent labelings to derive its ZDD. As described above, a nonequivalent labeling for a given supercell with the sublattice L M is defined as the maximum labeling in the lexicographical order among its equivalent labelings for the permutation group Σ M . In other words, a nonequivalent labeling is a labeling larger than any of its permuted structures for the permutation group. For permutation σ ∈ Σ M , the set of larger or equal labelings in the lexicographical order is expressed as Therefore, the set of the nonequivalent labelings C M,2 is then given by the intersection of C (σ) M,2 for all permutations as The intersection of ZDDs is efficiently obtained in a top-down manner [39] during the construction of an isomorphism-eliminated ZDD, C M,2 , using a frontierbased method [26], as described in Appendix B.

Multicomponent system
We extend the procedure for deriving ZDDs to multicomponent systems, keeping the binary structure of ZDD. First, we introduce an encoding of labeling Labeling c and one-hot encodingc have a one-to-one correspondence. The sum of k elements for site i in the onehot encoding must be one, because any one of k atomic species occupies site i. This means that the one-hot encoding must satisfy the one-of-k constraint of We then reformulate nonequivalent labelings using the one-hot encoding. A nonequivalent one-hot encoding can be defined as the largest one-hot encoding among its equivalent one-hot encodings in the lexicographical order. Therefore, the set of nonequivalent one-hot encod-ingsC M,k is written as where the action of permutation σ on one-hot encoding c, σ(c), is expressed as From this definition of the nonequivalent one-hot encodings, a ZDD of nonequivalent one-hot encodings is the intersection of isomorphism-eliminated ZDDs and a oneof-k ZDD representing whether the constraint is satisfied or not. Isomorphism-eliminated ZDD for permutation σ, C (σ) M,k , and one-of-k ZDDC one-of-k are written as and respectively. One-of-k ZDDC one-of-k as illustrated in Fig. 5 (a) is easily derived. Finally, the set of nonequivalent one-hot encodings is given as Note that if we impose additional constraints indicating prior knowledge on derivative structures such as energetically prohibited structures, a ZDD satisfying the additional constraints is derived by the intersection of the ZDDC M,k and the additional constraint ZDDs. Figure 5 (b) shows a ZDD representing ternary nonequivalent one-hot encodings of the two-dimensional supercell that is the same as that used in Section III C. Pólya's counting theorem or Eqn. (13) indicates that there exists 21 nonequivalent labelings. Because a 1path represents a nonequivalent one-hot encoding, the ZDD has 21 1-paths. The relationship between paths and nonequivalent one-hot encoding can be seen in the following example. Labeling c = (2, 1, 0, 0) is encoded to a one-hot representation as c = (c 4,2 ,c 4,1 ,c 4,0 ,c 3,2 ,c 3,1 ,c 3,0 ,c 2,2 ,c 2,1 ,c 2,0 ,c 1,2 ,c 1,1 ,c 1,0 ) = (1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1), (25) which is stored in the ZDD as the following 1-path. Note that the ZDD shown in Fig. 5 (b) is much more compact than the corresponding binary decision tree. The width of the ZDD or the maximum number of nodes corresponding to the same variable is only three, which corresponds to the number of nodes forc 3,0 ,c 2,2 , andc 2,1 .
On the other hand, the width of the binary decision tree is k |D M | = 3 4 = 81.
respectively. Then, we compare the performance of the present ZDD-based method with that of the previous method [15,16,18] implemented in enumlib [43]. As performed in the present ZDD-based method, the previous method (enumlib) enumerates nonequivalent labelings only for non-isomorphic permutation groups for a given index. We enumerate nonequivalent labelings to the limit of the index due to the computational resource [44]. Figure 8 shows the number of nonequivalent labelings for each non-isomorphic permutation group in (a) binary, (b) ternary, and (c) quaternary systems. As described before, the number of nonequivalent labelings is easily calculated by Pólya's counting theorem [12,13], and the number of nonequivalent labelings obtained from each ZDD coincides exactly with the number obtained by Pólya's counting theorem. The multiple numbers of nonequivalent labelings are found at most of the indexes in Fig. 8 because the number of nonequivalent labelings depends on the permutation group. Note that all nonequivalent labelings themselves can be obtained by tracing all paths in a ZDD, although we show only the number of nonequivalent labelings in Fig. 8.
For fcc, the present ZDD-based method enumerates nonequivalent labelings with up to 48, 31, and 26 sites for binary, ternary, and quaternary systems, respectively. The total number of derivative structures, which is the sum of the numbers of nonequivalent labelings over nonequivalent supercells, reaches approximately 10 17 in each of the binary, ternary, and quaternary systems. On the other hand, the total number of derivative structures is approximately 10 10 in every system when using the previous method implemented in enumlib, which is much smaller than that of the present ZDD-based method. For hcp, the number of sites in a derivative structure is twice its index. As can be seen in Fig. 8, the present ZDDbased method enumerates nonequivalent labelings with up to 50, 30, and 26 sites for the binary, ternary, and quaternary systems, respectively. The total number of derivative structures ranges approximately from 10 16 to 10 18 . Although the number of enumerated derivative structures can be too enormous to use them as candidates for optimization, it is expected to reduce candidates by introducing some useful constraints. Figure 9 shows the number of non-terminal nodes and the number of nonequivalent labelings in ZDDs. The diagonal line indicates that the number of non-terminal nodes is equal to that of nonequivalent labelings, which corresponds to a simple case of each labeling being expressed by a single node. Therefore, their ratio can be a simple estimation of the efficiency of a ZDD. As can be seen in Fig. 9, the number of non-terminal nodes is much smaller than that of nonequivalent labelings, which indicates that ZDD represents nonequivalent labelings efficiently. For example, ZDD compresses as many as approximately 10 12 nonequivalent labelings into approximately 10 8 non-terminal nodes.  Figure 10 shows the computational time required to enumerate derivative structures for a given index using the present ZDD-based method and the previous method (enumlib). The computational time required to construct a ZDD for each non-isomorphic permutation group is also shown in Fig. 10. In the ZDD-based method, the computational time required to enumerate derivative structures for a given index is the sum of the computational times required to construct ZDDs for non-isomorphic permutation groups. Similarly to the present method, the previous method (enumlib) enumerates nonequivalent labelings only for non-isomorphic permutation groups for a given index.
As can be seen in Fig. 10, both methods require exponential time with respect to the index. However, the two series of computational time indicate that the base of exponential time in the ZDD-based method is half that in the previous method. For example, the computational time is approximately 10 4 and 10 1 seconds at the index of 29 in the previous method and the ZDD-based method, respectively. The difference in the computational time at a larger index between the previous method and the ZDD-based method is expected to be much larger. Thus, the ZDD-based method to enumerate derivative structures is much more efficient than the previous method. In practice, memory consumption is also an essential aspect in enumerating derivative structures. A comparison of memory consumption between the ZDD-based method and the previous method is given in Appendix C.
Note that it is impossible to compare the present ZDD-based method and the previous method in a rigorous manner, because they have some differences in their procedures and implementations. However, the present discussion on their differences in terms of computational time order should remain valid. The main difference between their procedures is that the previous method (enumlib) excludes superperiodic structures and incomplete structures, but the present method includes such structures. Regarding the implementation of the methods, the present method and the previous method (enumlib) are implemented in C++ and Fortran, respectively.

VI. CONCLUSION
We have proposed an efficient procedure with a compact data structure of ZDD to enumerate derivative structures or nonequivalent labelings for given lattice and sites. We have applied the ZDD-based procedure to the enumeration of binary, ternary, and quaternary derivative structures from the simple fcc and hcp structures. The present approach significantly increases the possible number of derivative structures to be enumerated, which is as many as approximately 10 14 . Although the number of enumerated derivative structures can be too huge to perform structure optimization from the derivative structures, such a structure optimization will be possible by introducing ZDDs representing some useful constraints. The present ZDD-based method can be easily applied to the structure enumeration derived from the other structures. In addition to the derivative structure enumeration, ZDD and similar approaches should be powerful tools for solving combinatorial problems in physics and materials science that can be reformulated as the enumeration of subgraphs from a given graph. M,2 , for permutations π1 and π2. The solid and dotted arrows indicate 1-edges and 0-edges, respectively. The order of variables is fixed as c3, c2, and c1. To show the development process of ZDDs, partially determined labeling and its permuted labeling are also shown for each edge. The upper and lower blocks indicate (c3, c2, c1) and (c π(3) , c π(2) , c π(1) ), respectively. The asterisks indicate discarded labels that are not needed in a later comparison process.
terminal node has the value of one. An example of such a branching is c = (1 0 •) and π(c) = (0 • 1) for permutation π 1 . (3) If neither of the above relationships is decided, the labels used in the comparison are removed from the set of frontier labels that should continue to be compared. They do not need to be compared in a later process after comparing label c i and permuted label c π(i) for site i. Its example is c = (0 0 •) and π(c) = (0 • 0) for permutation π 1 . Labels c 1 and c π(1) are then discarded, and we denote them as c = ( * 0 •) and π(c) = ( * • 0), where * indicates a discarded label. (4) Two nodes are merged when the values of their frontier labels coincide [46].