BFS-Based Canonical Codes for Generating Graphs with Constraint Programming

In combinatorics and graph theory, graph enumeration involves systematically searching for all graphs that satisfy some given properties (such as a specific number of edges, vertex degree constraints, or structural features like girth). This has numerous applications in fields like chemical informatics and drug discovery, where molecules are modeled as graphs [5, 19]. To prevent us from generating isomorphic graphs, i.e., graphs that are equivalent up to a renaming of their vertices, we may use canonical codes, i.e., words associated with graphs such that two graphs are isomorphic if and only if they have the same canonical code: instead of enumerating graphs, we enumerate canonical codes, thus ensuring that all generated graphs are non-isomorphic.

Canonical codes may be based on adjacency matrices, i.e., 2-dimensional arrays $M$ such that $M_{ij}=1$ if $(i,j)$ is an edge, and 0 otherwise. A code may be obtained from $M$ by concatenating its rows (or columns), thus obtaining a binary word. Since different vertex permutations may produce different adjacency matrices and, therefore, different codes, the lexicographically smallest one is selected as the canonical code, thus uniquely identifying the graph. To avoid exhaustively considering all vertex permutations, symmetry-breaking constraints may be used. In [4], Codish et al. introduce symmetry-breaking constraints that enforce a non-descending row order in the adjacency matrix and ensure minimality under vertex permutations within specific partitioned sets.

The Nauty algorithm [13] generates a canonical labeling through an iterative partition refinement process, where vertices are colored based on their connectivity. Once a discrete partition is reached, it serves as the canonical label. Although there is no direct work that enumerates all canonical codes using Nauty’s labeling as constraints, as far as we know, the concept of symmetry breaking via structural graph information has been applied by Codish et al. [3]. They encode partition refinement into constraints, and introduce additional constraints to maintain minimality in the adjacency matrix while preserving partition structure. This approach provides insights into constructing compact and complete symmetry-breaking constraints, effectively reducing the need to evaluate all vertex permutations. Further studies on complete symmetry breaking are explored in [8, 9].

In addition to the adjacency matrix, Codish et al. demonstrate in [10] how various higher-dimensional graph invariants, such as vertex degrees and the cardinality of common neighbors, can be leveraged to define symmetry-breaking constraints. They compare different combinations of these invariants to evaluate their effectiveness in reducing the search space.

Beyond static symmetry breaking, which relies on a predefined set of permutations, Kirchweger and Szeider introduce in [11] a novel SMT-based approach for graph generation using dynamic symmetry breaking. The idea is to detect symmetries in the partially generated graph to enforce adjacency matrix minimality during the solving process, thus significantly improving efficiency compared to static symmetry-breaking methods.

Another class of canonical codes relies on graph traversals, where codes are sequences of traversed edges. Graph mining algorithms, such as Gaston [15] and gSpan [20], are usually based on Depth-First-Search (DFS). In [18] and [16], canonical codes based on Breadth-First-Search (BFS) are used for particular classes of graphs (Deterministic Finite Automata in [18] and hexagonal graphs in [16]) for which graph isomorphism may be solved in polynomial-time. In [14], some BFS-based symmetry-breaking predicates are introduced for generating connected graphs.

In this paper, we extend the work of [14] and introduce new symmetry breaking constraints for BFS-based canonical codes. Indeed, BFSs naturally construct spanning trees, and the orbits formed by the vertices within a spanning tree can be exploited to derive compact symmetry breaking constraints, that are efficiently handled by CP solvers. Also, every prefix of our canonical code is a canonical code, thus enabling the construction of dynamic symmetry-breaking constraints, similar to the approach in [11].

The paper is organized as follows. In Section 2, we introduce notations. In Section 3, we introduce BFS-based codes, and define a canonical BFS-based code as the smallest possible code. In Section 4, we describe a CP model for generating BFS-based codes. In Sections 5 and 6, we study properties of canonical codes. In Section 7, we introduce a global constraint that exploits these properties to ensure canonicity. In Section 8, we evaluate our method for generating connected claw-free cubic graphs, and compare it with the approach of [10], which relies on an adjacency matrix representation, as well as with the SMT-based approach of [11] and with Nauty [13]. We also introduce a more compact graph representation for this benchmark, thus enabling us to solve larger instances.

We note $[i,j]$ the set of all integers ranging from $i$ to $j$, $\mathrm{\#}S$ the cardinality of a set $S$, and ${⪯}_{lex}$ the lexicographic order for comparing sequences of integer values.

Two graphs $G_1=(V_1,E_1)$ and $G_2=(V_2,E_2)$ are isomorphic if there exists a bijection $f:V_1\to V_2$ that preserves edges, i.e., $\forall u,v\in V_1,(u,v)\in E_1\iff (f(u),f(v))\in E_2$. This problem belongs to $𝒩𝒫$, but it is not known to be in $𝒫$ nor to be $𝒩𝒫$-complete, and it is conjectured to be $𝒩𝒫$-intermediate.

Throughout this paper, we consider a connected graph $G=(V,E)$ such that $\mathrm{\#}V=n+1$ and $\mathrm{\#}E=n+m$, so that $0\le m\le \frac{(n-2)(n-1)}{2}$. We note $d_i$ the degree of a vertex $i\in V$.

A code is a sequence of $n+2m$ integer values, generated through a BFS of $G$ as described in Algorithm 1.

We use a FIFO queue $q$ to store the vertices that have been discovered but not yet treated, and we build two sequences $c_f$ and $c_b$ of $n$ and $2m$ integer values, respectively. When a vertex $v$ is pushed in $q$, we store in $num[v]$ its $q$-number, which ranges from 0 for the first pushed vertex to $n$ for the last pushed one. At each iteration of the while loop (Lines 4-13), a vertex $u$ is popped from $q$, and we consider each vertex $v$ adjacent to $u$.

• $\mathrm{■}$

  If $v$ is reached for the first time (Lines 7-8), $(num[u],num[v])$ is said to be a forward edge. $num[u]$ is the parent of $num[v]$, and it is added at the end of $c_f$. At the end of Algorithm 1, there are $n$ forward edges which form a spanning tree of $G$.

• $\mathrm{■}$

  If $v$ has already been reached (Lines 9-10), $(num[u],num[v])$ is said to be a Backward edge, and $num[u]$ and $num[v]$ are added at the end of $c_b$. To avoid treating twice a same edge, we add the condition that $num[u]$ must be smaller than $num[v]$.

The final code is obtained by concatenating $c_f$ and $c_b$ (Line 14). Hence, a code is a sequence of $n+2m$ integer values $p_1,p_2,\mathrm{\dots },p_n,x_1,y_1,x_2,y_2,\mathrm{\dots }{x}_m,y_m$ such that for each $i\in [1,n]$, $(p_i,i)$ is a forward edge and for each $i\in [1,m]$, $(x_i,y_i)$ is a backward edge. The first part of a code, from $p_1$ to $p_n$ is called forward code, whereas the second part of the code, from $x_1$ to $y_m$ is called backward code. When there is no ambiguity, commas between vertices in codes may be removed.

Different codes may be built for a same graph, depending on (i) which vertex is chosen first (Line 2) and (ii) the order considered to visit successors of $u$ (Line 6). In order to make the algorithm deterministic, we add a parameter $\pi$ which is a permutation of $V$ used to break ties: Line 2, the chosen vertex $v_0$ is the first vertex of $\pi$; Line 6, successors of $u$ are visited in the order defined by $\pi$. The code obtained when using $\pi$ to break ties is denoted BFS${}_G{}^{}(\pi )$.

Given a code, we can build the corresponding graph as the code allows us to reconstitute all its edges. For example, given BFS${}_G{}^{}(eabcd)=0001,12,13,24,34$, we reconstitute the set of edges $\{(0,1),(0,2),(0,3),(1,4),(1,2),(1,3),(2,4),(3,4)\}$.

As shown in Example 1, there exist different possible codes for $G$, depending on the considered permutation $\pi$. We define a total order on the set of all possible codes that may be associated with $G$ by considering a lexicographic order. Among all the possible codes for $G$, the smallest one according to this order is called the canonical code of $G$ and it is unique.

In the example of Figure 1, the canonical code is $0000,12,13,24,34$.

An important property to allow an efficient enumeration of canonical codes is that any prefix of a canonical code is a canonical code: this property allows us to stop completing a code whenever its prefix is not canonical.

We propose to generate codes with CP, thus allowing one to easily add other application-dependent constraints. Our CP model has the following integer variables:

• $\mathrm{■}$

  for each $i\in [1,n]$, $p_i$ corresponds to the parent of vertex $i$ in the spanning tree (in other words, $(p_i,i)$ is a forward edge), and its domain is $[0,i-1]$ because the parent of $i$ must have been discovered before $i$;

• $\mathrm{■}$

  for each $j\in [1,m]$, $x_j$ and $y_j$ correspond to the two endpoints of the $j$th backward edge and their domain is $[1,n]$ because an edge incident to $0$ cannot be backward;

• $\mathrm{■}$

  for each $i\in [0,n]$, $d_i$ corresponds to the degree of vertex $i$, and its domain is $[1,n-1]$.

The constraints are listed in Fig. 2. Constraints C1 to C7 are satisfied by any BFS-based code, even if it is not canonical:

• $\mathrm{■}$

  C1 is a consequence of the fact that the vertex whose $q$-number is 1 has been pushed in $q$ just after 0 and, therefore, its parent can only be 0.

• $\mathrm{■}$

  C2 and C4 are consequences of the fact that, at each iteration of the while loop (Lines 4-13), the vertex $u$ popped from $q$ has a $q$-number increased by one (as vertices are pushed in $q$ by increasing $q$-number).

• $\mathrm{■}$

  C3 is a consequence of the condition Line 9.

• $\mathrm{■}$

  C5 expresses the fact that, for each backward edge $(x_i,y_i)$, the parent of $y_i$ has been discovered before $x_i$. Indeed, let us suppose that this is not the case, i.e., $p_{y_i}>x_i$. In this case, the parent of $y_i$ would necessarily be $x_i$ because $y_i$ would not yet have been discovered when exploring the vertices adjacent to $x_i$ (Lines 7-8).

• $\mathrm{■}$

  C6 and C7 relate degrees with edge variables.

Constraints C8 to C10 prevent us from generating some non-canonical codes, i.e., codes that are not the smallest possible ones:

• $\mathrm{■}$

  C8 comes from the fact that $c$ starts with $d_0$ occurrences of $0$: if there exists a vertex $i$ such that $d_i>d_0$, then a smaller code is obtained by starting BFS from the vertex.

• $\mathrm{■}$

  If C9 is not satisfied, then a smaller code is obtained by exchanging $x_i{y}_i$ and $x_{i+1}{y}_{i+1}$.

• $\mathrm{■}$

  if C10 is not satisfied, then a smaller code is obtained by visiting $i+1$ before $i$.

We have experimentally compared this first CP model with the model introduced in [14] on a toy problem that aims at enumerating all graphs with $k$ vertices and $2k-2$ edges. Both models have been implemented in Choco, and compute the same sets of solutions. However, our model is more efficient: when $k=5$ (resp. 6, 7, and 8), it needs 0.02 (resp. 0.1, 0.4, and 4.3) seconds whereas the model of [14] needs 0.04 (resp. 0.2, 1.4, and 28.4) seconds. This comes from the fact that (i) we order sibling vertices with respect to the number of children (thanks to constraint C10) instead of sub-tree weights, and (ii) we explicitly model backward edges with $x_j$ and $y_j$ variables instead of using an adjacency boolean matrix.

Constraints C1 to C10 allow us to enumerate all canonical BFS-based codes, but some non-canonical codes may also be enumerated. Hence, in the next two sections we introduce two additional properties that must be satisfied by canonical codes, and that are used to propagate a global constraint that ensures the canonicity of BFS-based codes, as described in Section 7.

In this section, we introduce a property which is used in Section 7 to ensure the canonicity of prefixes of forward codes. Forward code prefixes correspond to trees, and the tree isomorphism problem may be solved in polynomial time by computing vertex labels [1]. We exploit similar labels but, unlike the original algorithm of [1], we do not rename labels with integer values at each level of the tree, but define the label of a vertex $i$ as a sequence that contains an integer value for each vertex in the subtree rooted in $i$. This allows us to compute labels in a more incremental way during the search.

For example, let us consider the tree $T$ displayed on the left of Fig. 3. The forward code associated with this tree is $c_f=000113$. When $i=1$, we have ${\mathrm{\Gamma }}_1=145$ and $l_1=200$ because $\mathrm{\#}{\mathrm{𝐶ℎ}}_{c_f}(\mathit{1})=\mathit{2}$ and $\mathrm{\#}{\mathrm{𝐶ℎ}}_{c_f}(\mathit{4})=\mathrm{\#}{\mathrm{𝐶ℎ}}_{c_f}(\mathit{5})=\mathit{0}$. When $i=0$, we have ${\mathrm{\Gamma }}_0=0123456$ and $l_0=3201000$ because $\mathrm{\#}{\mathrm{𝐶ℎ}}_{c_f}(\mathit{0})=\mathit{3}$, $\mathrm{\#}{\mathrm{𝐶ℎ}}_{c_f}(\mathit{1})=\mathit{2}$, $\mathrm{\#}{\mathrm{𝐶ℎ}}_{c_f}(\mathit{3})=\mathit{1}$, and $\mathrm{\#}{\mathrm{𝐶ℎ}}_{c_f}(\mathit{2})=\mathrm{\#}{\mathrm{𝐶ℎ}}_{c_f}(\mathit{4})=\mathrm{\#}{\mathrm{𝐶ℎ}}_{c_f}(\mathit{5})=\mathrm{\#}{\mathrm{𝐶ℎ}}_{c_f}(\mathit{6})=\mathit{0}$.

Given the label $l_0$ of the root of a tree, we can reconstitute this tree as each value at position $i$ in $l_0$ gives the number of children of vertex $i$. Hence, there is a bijection between each forward code and the label of the root of its associated tree. For example, we display in Fig. 3 two isomorphic trees $T$ and $T^{\prime }$. $T$ corresponds to the forward code $000113$ and the label of its root is $3201000$, whereas $T^{\prime }$ corresponds to the forward code $000112$ and the label of its root is $3210000$.

An interesting property of vertex labels is that the children of a vertex in the tree associated with a canonical forward code prefix have non-decreasing labels, as stated below.

Given a forward code prefix, we propose to exploit automorphisms in the tree associated with this prefix in order to break symmetries on backward edges. An automorphism of a tree $T=(V_T,E_T)$ is a permutation $\sigma$ of $V_T$ which is an isomorphism from $T$ to $T$. Given a forward code prefix $c_f$, we only consider automorphisms that are still automorphisms on all possible extensions of $c_f$ to full canonical forward codes, as defined below.

For example, let us consider the tree $T^{\prime }$ displayed in Fig. 3. The canonical code of $T^{\prime }$ is $000112$. Let us assume that $n=9$, i.e., we still have to add three forward edges to the tree. Let ${\sigma }_{45}$ be the automorphism that only exchanges vertices 4 and 5 and leaves unchanged all other vertices, i.e., ${\sigma }_{45}(4)=5,{\sigma }_{45}(5)=4$, and $\forall i\notin \{4,5\},{\sigma }_{45}(i)=i$. ${\sigma }_{45}$ does not belong to $Aut(000112)$. Indeed, there exist extensions of $000112$ for which ${\sigma }_{45}$ is not an automorphism. For example, if we add the forward edges $(3,7)$ and $(4,8)$ to $T^{\prime }$, leading to the forward code $00011234$, then ${\sigma }_{45}$ is no longer an automorphism because the subtree rooted in 4 is no longer isomorphic to the subtree rooted in 5.

Now, let us consider the tree obtained from $T^{\prime }$ by adding the forward edge $(6,7)$. The associated forward code is $0001126$. In this case, ${\sigma }_{45}$ belongs to $Aut(0001126)$ because it is still an automorphism on all canonical extensions of $0001126$. Indeed, 4 and 5 cannot be parent (because 6, which is at the same level as 4 and 5 but with a larger q-number, already has a child).

The following property allows us to exploit automorphisms to break symmetries on backward edges.

Properties 1 and 2 are used to propagate the global constraint defined below.

Cubic graphs are graphs in which each vertex has degree 3. A cubic graph is claw-free if it does not contain $K_{1,3}$ as an induced subgraph. This is equivalent to requiring that every vertex must participate in at least one triangle.

To enumerate all claw-free cubic graphs, we define a CP model composed of Constraint C1 to C10 combined with our global cc constraint. To enforce the cubic degree condition, we set the domain of each degree variable $d_i$ to $\{3\}$, for each $i\in [0,n]$.

To ensure claw-freeness, we introduce a global constraint that ensures that each vertex is involved in at least one triangle. The propagator of this constraint maintains two sets of edges as proposed in [2]: mandatory edges, which must be in the solution, and possible edges, which may be included. We use sparse sets to efficiently maintain these sets [12].

A dynamic variable selection strategy is employed to ensure that $de{p}_p-de{p}_x>1$ is maintained throughout the search, where $de{p}_p$ and $de{p}_x$ are the depth levels of the last instantiated forward and backward edge variables, respectively. By prioritizing the instantiation of $p$-variables while interleaving backward edge variables, we avoid the enumeration of useless spanning trees. Additionally, for every $i\in [1,m]$, $y_i$ is instantiated just after $x_i$.

As the procedure for checking the canonicity is rather expensive, we introduce a parameter $f$ which allows us to control the frequency of canonicity checking: When $f=1$, canonicity is checked after each edge assignment; when $f>1$, it is checked every $f$ edge assignment. Of course, when all variables are assigned, the canonicity check is performed, whatever the value of $f$ is, in order to ensure that the global constraint is satisfied. A similar parameter is used in the SMT-based approach of [11].

Table 2 shows the results for $n\in [20,44]$ by steps of 2, when our parameter $f$ belongs to $\{1,10,20,30\}$. When $f=1$, run times are very often longer than when $f\ge 10$. A good tradeoff is reached when $f\in \{10,20\}$.

In Table 2, we also display the results reported by Codish et al. in [10] (the code of this approach is not available). This approach does not break all symmetries and, therefore, it may compute redundant graphs that are isomorphic to previously enumerated graphs. For example, when $n=32$, there are $731$ different graphs whereas the approach of [10] computes $29,069$ solutions. As a consequence, this approach does not scale well and cannot be used to solve larger instances within a reasonable amount of time. As a comparison, our method which only computes non-isomorphic graphs achieves a speed-up of more than 60 for $n\ge 30$ (note however that the two approaches have been run on different computers).

We have adapted the SMT-based approach described in [11] to generate claw-free cubic graphs, and we display the results obtained with this approach, on the same computer as the one used in our experiments. To handle the claw-free constraint, we directly utilize the solver’s built-in feature, "–forbidden-induced-subgraphs", to prohibit the predefined induced subgraph $K_{1,3}$. In [11], a parameter similar to $f$ is used to control the frequency of "minimality check", a non-polynomial operation that verifies whether the partially defined graph is canonical in their approach. We report results obtained with $f\in \{1,10,20,30\}$ (the default value for this parameter is 20). The best results of SMT are obtained when $f=10$, and our method consistently outperforms SMT across all instances.

Finally, we display results obtained with Nauty [13] (using geng -F to ensure claw-freeness). Our approach is always more efficient than Nauty, but it needs more memory (e.g., 444584Kb instead of 22760Kb when $n=30$).

We introduce a new canonical graph encoding based on a BFS traversal, which is well-suited for constraint-based approaches. Using this encoding, we define a CP model that enumerates all non-isomorphic graphs satisfying specific constraints by generating their corresponding canonical codes. The key aspect of this encoding is identifying a spanning tree that yields the lexicographically smallest representation of the graph. We demonstrate how this encoding enables the formulation of fundamental static symmetry-breaking constraints. Additionally, we leverage the crucial property that every prefix of a canonical code must also be canonical, allowing us to define a global constraint that dynamically enforces symmetry-breaking, ensuring the canonicity of the generated codes.

We tested our approach on connected claw-free cubic graphs, which are regular and inherently challenging for BFS-based methods. However, our experimental results demonstrate that our approach outperforms existing state-of-the-art techniques based on adjacency matrix representations. Beyond this, we explored the properties of local structures within this class of graphs that allow us to contract graphs into more compact labeled graphs. Our BFS-based canonical code may be easily extended to labeled graphs, and this allows us to enumerate graphs of larger order, up to 60 vertices. As part of our future work, we aim to refine our encoding for such labeled graphs and develop dedicated propagation strategies.

Our approach is well-suited for enumerating graphs with a limited diameter, such as diameter-2-critical graphs and required-girth extremal graphs [11]: In these cases, the number of canonical spanning trees remains manageable, allowing them to be precomputed and used for static symmetry breaking. We plan to publish further results on these aspects in the near future, as they could not be included in this work due to space limitations.