pipeline

def compile_program(
    prepared: SamplingGraph,
    *,
    mode: DecompositionMode,
    strategy: DecompositionStrategy = "cat5",
) -> CompiledProgram:
    """Compile a prepared graph into an executable sampling program.

    This function performs the second phase of compilation:
    1. Split the graph into connected components
    2. For each component:
       - Plug outputs according to mode (sequential or joint)
       - Reduce each plugged graph
       - Perform stabilizer rank decomposition
       - Compile into CompiledScalarGraphs objects
    3. Assemble into CompiledProgram with output ordering

    Args:
        prepared: The prepared graph from prepare_graph().
        mode: Decomposition mode:
            - "sequential": For sampling - creates [0, 1, 2, ..., n] circuits
            - "joint": For probability estimation - creates [0, n] circuits
        strategy: Stabilizer rank decomposition strategy.
            Must be one of "cat5", "bss", "cutting".

    Returns:
        A CompiledProgram ready for sampling.

    """
    components = connected_components(prepared.graph)

    # Determine global f-indices (numerically sorted) from the prepared graph
    f_indices_global = _get_f_indices(prepared.graph)
    num_outputs = prepared.num_outputs

    direct_entries: list[tuple[int, int, bool]] = []  # (output_idx, f_idx, flip)
    compiled_components: list[CompiledComponent] = []
    compiled_output_order: list[int] = []

    sorted_components = sorted(components, key=lambda c: len(c.output_indices))

    for component in sorted_components:
        direct = classify_direct(component)
        if direct is not None:
            f_idx, flip = direct
            direct_entries.append((component.output_indices[0], f_idx, flip))
        else:
            compiled = _compile_component(
                component=component,
                f_indices_global=f_indices_global,
                mode=mode,
                strategy=strategy,
            )
            compiled_components.append(compiled)
            compiled_output_order.extend(component.output_indices)

    # Sort direct entries by output index so that — together with the output
    # prioritisation in transform_error_basis — the concatenated layout often
    # matches the original output order, sparing a reindex at sample time.
    direct_entries.sort()
    direct_output_order = [e[0] for e in direct_entries]
    direct_f_indices = [e[1] for e in direct_entries]
    direct_flips = [e[2] for e in direct_entries]

    output_order = np.array(direct_output_order + compiled_output_order, dtype=np.int32)
    reindex = np.argsort(output_order)
    is_identity = np.array_equal(reindex, np.arange(len(output_order)))

    return CompiledProgram(
        components=tuple(compiled_components),
        direct_f_indices=jnp.array(direct_f_indices, dtype=jnp.int32),
        direct_flips=jnp.array(direct_flips, dtype=jnp.bool_),
        output_order=jnp.asarray(output_order),
        output_reindex=None if is_identity else jnp.asarray(reindex),
        num_outputs=num_outputs,
        num_detectors=prepared.num_detectors,
    )