Extract minimum assembly pathways as directed acyclic graphs of fragments

[jdaymude]

Resolves #22, enabling minimum assembly pathway extraction. Subsumes #146 and #147.

There's a lot to talk about here, but the TL;DR is this: we now finally support assembly pathway reconstruction, we do so in a reasonably efficient way (5–20% performance overhead to search when reconstructing pathways, negligible performance differences if pathway reconstruction is disabled), and our pathway data format is more informative and interoperable with downstream tools than what's out there right now.

Specifying Assembly Pathway Reconstruction

There are several interface-level changes:

• the CLI now has a --pathways option
• the Rust index_search function now takes a mandatory max_pathways: Option<usize> parameter
• the Python index_search function now takes a max_pathways: int/None parameter (default None)

All of these basically work the same way:

• If --pathways is not set (CLI) or max_pathways = None (Rust/Python), assembly pathway reconstruction is disabled and things mostly work as before.
• If --pathways n (CLI) or max_pathways = n (Rust/Python) for some integer n > 0, we return at most n minimum assembly pathways discovered during assembly index search.
• If --pathways 0 (CLI) or max_pathways = 0 (Rust/Python), we return all minimum assembly pathways discovered during assembly index search.

Note that the last of these settings does not necessarily return all possible minimum assembly pathways. When search optimizations like bounding and memoization are enabled, we prune any search states that cannot yield a better bound on the assembly index. This includes assembly states that may attain equally good assembly index bounds, and in particular, those that would attain the true assembly index after one such state is already found. These states may also potentially yield minimum assembly pathways, but they are pruned for search efficiency.

Assembly Pathway Representation

Assembly pathways are represented as directed, acyclic multigraphs (DAGs) whose nodes are molecular graph fragments and whose edges represent joining operations. Our data structure is similar to what you would need for the Croninlab visualization, but with a few additional improvements. Formally:

• Each DAG node represents a canonical fragment (i.e., a connected subgraph) of the target molecular graph. These fragments are "canonical" in the sense that they are deduplicated by isomorphism. For example, if an assembly pathway joins two CC single bonds together to form a CCC and then joins another single CC bond to that structure to form a CCCC chain, the CC fragment will only be represented once as a node in the DAG.
• The DAG contains directed edges $u_1 \to v$ and $u_2 \to v$ if and only if some fragment $f_1$ isomorphic to $u_1$ and some fragment $f_2$ isomorphic to $u_2$ are joined to produce a fragment $f_3$ isomorphic to $v$ in the assembly pathway. This means that an edge $u \to v$ may appear in the DAG multiple times (hence it is a multigraph) if, e.g., two isomorphic fragments are joined to each other.
• Digging deeper, each DAG node and edge has a label that specifies exactly which bonds from the target molecular graph its fragment contains. For DAG nodes, this label specifies the bonds of the first real instance of that node's canonical fragment. For DAG edges $u_1 \to v$ and $u_2 \to v$ representing the joining of fragment $f_1$ isomorphic to $u_1$ and fragment $f_2$ isomorphic to $u_2$, the edges are labeled with the bonds of $f_1$ and $f_2$, respectively. This means that the pathway not only represents which canonical fragments are joined to produce other canonical fragments, but which specific instances of those canonical fragments were used in joining operations.

Pathways are returned as follows:

• The CLI prints the target molecular graph and any reconstructed assembly pathways as DOT-formatted strings. For the pathway graphs, node and edge labels are sets of bond indices corresponding to the bonds in the target molecular graph.

Example Output for Anthracene

> ./target/release/assembly-theory data/checks/anthracene.mol --parallel none --pathways 0
Assembly Index:  6
Matches:         466
States Searched: 491
Pathways Found:  2

Molecule Graph: graph {
    0 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    1 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    2 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    3 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    4 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    5 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    6 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    7 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    8 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    9 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    10 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    11 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    12 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    13 [ label = "Atom { element: Carbon, capacity: 0 }" ]
    8 -- 7 [ label = "Double" ]
    7 -- 4 [ label = "Single" ]
    4 -- 3 [ label = "Double" ]
    3 -- 0 [ label = "Single" ]
    0 -- 1 [ label = "Double" ]
    1 -- 2 [ label = "Single" ]
    2 -- 5 [ label = "Double" ]
    8 -- 9 [ label = "Single" ]
    9 -- 10 [ label = "Single" ]
    10 -- 11 [ label = "Double" ]
    11 -- 12 [ label = "Single" ]
    12 -- 13 [ label = "Double" ]
    13 -- 8 [ label = "Single" ]
    4 -- 5 [ label = "Single" ]
    5 -- 6 [ label = "Single" ]
    6 -- 9 [ label = "Double" ]
}

Pathway 0: digraph {
    0 [ label = "{14}" ]
    1 [ label = "{15}" ]
    2 [ label = "{14, 15}" ]
    3 [ label = "{7, 14, 15}" ]
    4 [ label = "{6, 7, 14, 15}" ]
    5 [ label = "{3, 4, 5, 6, 7, 14, 15}" ]
    6 [ label = "{2, 3, 4, 5, 6, 7, 13, 14, 15}" ]
    7 [ label = "{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}" ]
    0 -> 2 [ label = "{14}" ]
    1 -> 2 [ label = "{15}" ]
    0 -> 3 [ label = "{7}" ]
    2 -> 3 [ label = "{14, 15}" ]
    1 -> 4 [ label = "{6}" ]
    3 -> 4 [ label = "{7, 14, 15}" ]
    4 -> 5 [ label = "{6, 7, 14, 15}" ]
    3 -> 5 [ label = "{3, 4, 5}" ]
    2 -> 6 [ label = "{2, 13}" ]
    5 -> 6 [ label = "{3, 4, 5, 6, 7, 14, 15}" ]
    6 -> 7 [ label = "{2, 3, 4, 5, 6, 7, 13, 14, 15}" ]
    5 -> 7 [ label = "{0, 1, 8, 9, 10, 11, 12}" ]
}

Pathway 1: digraph {
    0 [ label = "{7}" ]
    1 [ label = "{15}" ]
    2 [ label = "{7, 15}" ]
    3 [ label = "{7, 14, 15}" ]
    4 [ label = "{6, 7, 14, 15}" ]
    5 [ label = "{3, 4, 5, 6, 7, 14, 15}" ]
    6 [ label = "{2, 3, 4, 5, 6, 7, 13, 14, 15}" ]
    7 [ label = "{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}" ]
    0 -> 2 [ label = "{7}" ]
    1 -> 2 [ label = "{15}" ]
    0 -> 3 [ label = "{14}" ]
    2 -> 3 [ label = "{7, 15}" ]
    1 -> 4 [ label = "{6}" ]
    3 -> 4 [ label = "{7, 14, 15}" ]
    4 -> 5 [ label = "{6, 7, 14, 15}" ]
    3 -> 5 [ label = "{3, 4, 5}" ]
    2 -> 6 [ label = "{2, 13}" ]
    5 -> 6 [ label = "{3, 4, 5, 6, 7, 14, 15}" ]
    6 -> 7 [ label = "{2, 3, 4, 5, 6, 7, 13, 14, 15}" ]
    5 -> 7 [ label = "{0, 1, 8, 9, 10, 11, 12}" ]
}

• The Rust index_search function returns a Vec<Pathway>, where the Pathway struct represents a single assembly pathway as a petgraph::Digraph whose nodes/edges are labeled by bit_set::BitSets indexing their fragments' bonds, following the specifications above.
• The Python index_search function returns a list[str] containing the DOT-formatted string representations of all reconstructed assembly pathways, appearing exactly as they do in the CLI output.

This data format is particularly nice for the Python interface, which can load these DOT strings up in graphviz, parse the node and edge labels as sets using Python's ast.literal_eval, and then throw those sets into RDKit's PathToSubmol to get the actual molecular fragments used in the assembly pathway and visualize them. This PR does not implement this functionality (it's big enough as is), but sets up all the necessary data structures.

Performance Impact

When assembly pathway reconstruction is disabled, search performance is not impacted:

Benchmark 06e6582 (this PR) vs. a992612 (current main) as a baseline

bench_bounds/gdb13_1201/no-bounds
                        time:   [65.494 ms 65.956 ms 66.446 ms]
                        change: [+1.4671% +2.4459% +3.5088%] (p = 0.00 < 0.05)
                        Performance has regressed.
bench_bounds/gdb13_1201/log
                        time:   [63.690 ms 64.226 ms 64.755 ms]
                        change: [−0.7080% +0.3538% +1.3847%] (p = 0.52 > 0.05)
                        No change in performance detected.
bench_bounds/gdb13_1201/int
                        time:   [60.384 ms 60.987 ms 61.538 ms]
                        change: [−0.2027% +1.1940% +2.4425%] (p = 0.08 > 0.05)
                        No change in performance detected.
bench_bounds/gdb13_1201/int-vec
                        time:   [56.909 ms 57.341 ms 57.788 ms]
                        change: [−0.3523% +0.9907% +2.4137%] (p = 0.17 > 0.05)
                        No change in performance detected.
bench_bounds/gdb13_1201/int-matchable
                        time:   [61.348 ms 61.755 ms 62.203 ms]
                        change: [+0.8296% +1.9401% +3.0818%] (p = 0.00 < 0.05)
                        Change within noise threshold.
bench_bounds/gdb17_200/no-bounds
                        time:   [208.20 ms 210.56 ms 212.99 ms]
                        change: [−0.9079% +0.4494% +1.8667%] (p = 0.54 > 0.05)
                        No change in performance detected.
bench_bounds/gdb17_200/log
                        time:   [148.47 ms 149.23 ms 149.93 ms]
                        change: [−0.9271% −0.2737% +0.4107%] (p = 0.45 > 0.05)
                        No change in performance detected.
bench_bounds/gdb17_200/int
                        time:   [51.486 ms 51.710 ms 51.955 ms]
                        change: [+0.5067% +1.1016% +1.7279%] (p = 0.00 < 0.05)
                        Change within noise threshold.
bench_bounds/gdb17_200/int-vec
                        time:   [50.365 ms 50.619 ms 50.895 ms]
                        change: [+1.2835% +1.9808% +2.7218%] (p = 0.00 < 0.05)
                        Performance has regressed.
bench_bounds/gdb17_200/int-matchable
                        time:   [40.506 ms 40.661 ms 40.821 ms]
                        change: [+1.2662% +1.9492% +2.6835%] (p = 0.00 < 0.05)
                        Performance has regressed.
bench_bounds/checks/no-bounds
                        time:   [164.24 ms 167.38 ms 170.49 ms]
                        change: [−0.5238% +2.1960% +4.9069%] (p = 0.11 > 0.05)
                        No change in performance detected.
bench_bounds/checks/log time:   [94.569 ms 95.478 ms 96.432 ms]
                        change: [−1.9663% −0.6575% +0.6738%] (p = 0.35 > 0.05)
                        No change in performance detected.
bench_bounds/checks/int time:   [10.945 ms 11.002 ms 11.060 ms]
                        change: [−0.3730% +0.5608% +1.6315%] (p = 0.32 > 0.05)
                        No change in performance detected.
bench_bounds/checks/int-vec
                        time:   [7.8191 ms 7.8653 ms 7.9188 ms]
                        change: [−0.5047% +0.4176% +1.2699%] (p = 0.39 > 0.05)
                        No change in performance detected.
bench_bounds/checks/int-matchable
                        time:   [6.2968 ms 6.3299 ms 6.3710 ms]
                        change: [−0.5332% +0.4232% +1.4514%] (p = 0.43 > 0.05)
                        No change in performance detected.
bench_bounds/coconut_55/no-bounds
                        time:   [1.8370 s 1.8948 s 1.9590 s]
                        change: [−2.6141% +1.1449% +5.6358%] (p = 0.58 > 0.05)
                        No change in performance detected.
bench_bounds/coconut_55/log
                        time:   [1.3173 s 1.3410 s 1.3647 s]
                        change: [−0.1705% +2.2820% +4.7587%] (p = 0.09 > 0.05)
                        No change in performance detected.
bench_bounds/coconut_55/int
                        time:   [129.00 ms 130.41 ms 131.73 ms]
                        change: [−1.1622% +0.3864% +1.9602%] (p = 0.62 > 0.05)
                        No change in performance detected.
bench_bounds/coconut_55/int-vec
                        time:   [114.17 ms 115.79 ms 117.40 ms]
                        change: [−1.7570% −0.0124% +1.7974%] (p = 0.98 > 0.05)
                        No change in performance detected.
bench_bounds/coconut_55/int-matchable
                        time:   [43.933 ms 44.095 ms 44.265 ms]
                        change: [−0.2354% +0.3520% +0.9733%] (p = 0.28 > 0.05)
                        No change in performance detected.

bench_memoize/gdb13_1201/no-memoize
                        time:   [46.247 ms 46.503 ms 46.768 ms]
                        change: [+1.8958% +2.7615% +3.5646%] (p = 0.00 < 0.05)
                        Performance has regressed.
bench_memoize/gdb13_1201/nauty-index
                        time:   [63.671 ms 64.097 ms 64.520 ms]
                        change: [+0.8332% +1.9288% +3.0107%] (p = 0.00 < 0.05)
                        Change within noise threshold.
bench_memoize/gdb13_1201/tree-nauty-index
                        time:   [60.756 ms 61.208 ms 61.666 ms]
                        change: [+1.2603% +2.3707% +3.3642%] (p = 0.00 < 0.05)
                        Performance has regressed.
bench_memoize/gdb17_200/no-memoize
                        time:   [25.407 ms 25.547 ms 25.760 ms]
                        change: [−1.9346% +2.7001% +7.6541%] (p = 0.29 > 0.05)
                        No change in performance detected.
bench_memoize/gdb17_200/nauty-index
                        time:   [43.316 ms 43.473 ms 43.608 ms]
                        change: [+1.0217% +1.7233% +2.3440%] (p = 0.00 < 0.05)
                        Performance has regressed.
bench_memoize/gdb17_200/tree-nauty-index
                        time:   [40.503 ms 40.733 ms 40.941 ms]
                        change: [+1.5477% +2.3846% +3.2019%] (p = 0.00 < 0.05)
                        Performance has regressed.
bench_memoize/checks/no-memoize
                        time:   [4.2806 ms 4.3040 ms 4.3341 ms]
                        change: [−1.0161% −0.1236% +0.6600%] (p = 0.79 > 0.05)
                        No change in performance detected.
bench_memoize/checks/nauty-index
                        time:   [7.1944 ms 7.2393 ms 7.2811 ms]
                        change: [−0.2995% +0.4269% +1.1962%] (p = 0.28 > 0.05)
                        No change in performance detected.
bench_memoize/checks/tree-nauty-index
                        time:   [6.2994 ms 6.3144 ms 6.3357 ms]
                        change: [−0.8565% +0.1030% +1.0101%] (p = 0.83 > 0.05)
                        No change in performance detected.
bench_memoize/coconut_55/no-memoize
                        time:   [34.316 ms 34.437 ms 34.567 ms]
                        change: [−1.7318% −0.0664% +1.4859%] (p = 0.94 > 0.05)
                        No change in performance detected.
bench_memoize/coconut_55/nauty-index
                        time:   [48.713 ms 48.955 ms 49.221 ms]
                        change: [−0.2698% +0.4273% +1.0910%] (p = 0.24 > 0.05)
                        No change in performance detected.
bench_memoize/coconut_55/tree-nauty-index
                        time:   [44.061 ms 44.199 ms 44.339 ms]
                        change: [+0.6638% +1.0976% +1.5606%] (p = 0.00 < 0.05)
                        Change within noise threshold.

A new benchmark, bench_removal_orders, compares search performance across three scenarios: (1) assembly pathway reconstruction disabled, (2) reconstructing one minimum assembly pathway discovered during search, and (3) reconstructing all minimum assembly pathways discovered during search. The raw data is below, but search incurs a 5–20% overhead when collecting match removal orders for later assembly pathway reconstruction. The overhead is smaller when fewer match removal orders are collected (e.g., for coconut_55, the overhead is 11.84% for one match removal order but 15.94% for all).

Benchmark 06e6582 (this PR) on search with assembly pathway reconstruction enabled

bench_removal_orders/gdb13_1201/none
                        time:   [60.436 ms 60.894 ms 61.388 ms]
bench_removal_orders/gdb13_1201/one
                        time:   [63.645 ms 64.122 ms 64.705 ms]
bench_removal_orders/gdb13_1201/all
                        time:   [63.722 ms 64.253 ms 64.818 ms]
bench_removal_orders/gdb17_200/none
                        time:   [40.628 ms 40.794 ms 40.977 ms]
bench_removal_orders/gdb17_200/one
                        time:   [47.618 ms 47.897 ms 48.192 ms]
bench_removal_orders/gdb17_200/all
                        time:   [48.229 ms 48.536 ms 48.858 ms]
bench_removal_orders/checks/none
                        time:   [6.3210 ms 6.3424 ms 6.3767 ms]
bench_removal_orders/checks/one
                        time:   [7.1797 ms 7.2185 ms 7.2626 ms]
bench_removal_orders/checks/all
                        time:   [7.4138 ms 7.4457 ms 7.4872 ms]
bench_removal_orders/coconut_55/none
                        time:   [43.985 ms 44.134 ms 44.286 ms]
bench_removal_orders/coconut_55/one
                        time:   [49.143 ms 49.359 ms 49.580 ms]
bench_removal_orders/coconut_55/all
                        time:   [51.036 ms 51.170 ms 51.310 ms]

Algorithm Details

Tracking match removal orders during search was something we were already doing; #149 and #150 handled the necessary setup to be used for assembly pathway reconstruction. Our reconstruction algorithm is an efficient implementation of the duplicate/remnant algorithm described in the Supporting Information of Seet et al. (2025).