Spinful Fermi-Hubbard dynamics - when mode ordering meets hardware connectivity

import copy
import math
from collections import defaultdict
from importlib.metadata import version
from typing import NamedTuple

import ffsim
import matplotlib.pyplot as plt
from matplotlib.colors import to_rgb
from matplotlib.patches import FancyArrowPatch
import numpy as np
import rustworkx as rx
import scipy.sparse.linalg as spla
from rustworkx.visualization import mpl_draw
from qiskit import QuantumCircuit
from qiskit.circuit import Parameter, QuantumRegister
from qiskit.circuit.library import (
    RXGate,
    RZZGate,
)
from qiskit.passmanager import MultiStagePassManager
from qiskit.quantum_info import (
    Operator,
    SparseObservable,
)
from qiskit.synthesis import LieTrotter
from qiskit.synthesis.evolution import EvolutionSynthesis
from qiskit.transpiler import CouplingMap, InstructionProperties
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

from qiskit_ibm_runtime.fake_provider import FakeFez, FakeNighthawk

# Importing this module registers the rule that folds out-of-range rzz angles into (0, pi/2] (used by `WrapAngles`).
import qiskit_ibm_runtime.utils.backend_converter  # noqa: F401

from qiskit_fermions.circuit import FermionicCircuit
from qiskit_fermions.circuit.library import Evolution
from qiskit_fermions.circuit.library.synthesis import FermionicLieTrotter
from qiskit_fermions.mappers import map_transfer_vertex_generators
from qiskit_fermions.mappers.library import (
    fermion_operator,
)
from qiskit_fermions.operators import TransferVertexOperator
from qiskit_fermions.transpiler import (
    FermionicCircuitToDAG,
    QuantumDAGToCircuit,
)
from qiskit_fermions.transpiler.passes import (
    CustomF2QLayout,
    F2QSynthesis,
    MapperFnEvolutionSynthesis,
)

for package in (
    "qiskit",
    "qiskit-fermions",
    "qiskit-ibm-runtime",
    "ffsim",
    "rustworkx",
):
    print(f"{package:20s} {version(package)}")

L_TINY = 2  # small enough to print every term of an operator
L_SMALL = 4  # the size of every explicit check and drawing below (the spinless guide uses four sites, too)
NUM_SEEDS = 20  # the number of seeds to try for randomized transpilation runs
GATE_SETS = (
    "cz basis",
    "fractional",
)  # the two gate sets of the routing studies

# Colors of the drawings: one meaning per color throughout the notebook. Spinless sites and the two spin species are shades of
# one blue (spin up dark, spin down light); hops are a dark/light pair of orange (east/west); ancillas are green, and red is
# kept for links that are not couplers.
SPIN_UP, SPIN_DOWN, SPINLESS, ANCILLA = (
    "#1b4f8a",
    "#9ecae1",
    "tab:blue",
    "tab:green",
)
HOP_EAST, HOP_WEST = (
    "#e6550d",
    "#fdae6b",
)  # hops (transfer operators) along a chain, dark / light orange: pointing to the larger / the smaller mode index

# Device-specific coupling map constants
ROWS, WIDTH = 12, 10

# The style of every curve in the plots below, keyed by its label. Routed variants: solid = cz basis, dashed = fractional gates.
STYLES = {
    "LieTrotter, line": {"color": "0.7", "marker": "o"},
    "flow set, line": {"color": "0.3", "marker": "s"},
    "interleaved, 1 ancilla, Heron, cz basis": {
        "color": "tab:green",
        "marker": "o",
    },
    "interleaved, 1 ancilla, Heron, fractional": {
        "color": "tab:green",
        "marker": "o",
        "linestyle": "--",
    },
    "interleaved, 1 ancilla, Nighthawk, cz basis": {
        "color": "tab:olive",
        "marker": "^",
    },
    "interleaved, 1 ancilla, Nighthawk, fractional": {
        "color": "tab:olive",
        "marker": "^",
        "linestyle": "--",
    },
    "blocked, 2 ancillas, Nighthawk, cz basis": {
        "color": "tab:blue",
        "marker": "D",
    },
    "blocked, 2 ancillas, Nighthawk, fractional": {
        "color": "tab:blue",
        "marker": "D",
        "linestyle": "--",
    },
    "blocked, 2 ancillas, Heron, cz basis": {
        "color": "tab:red",
        "marker": "v",
    },
    "blocked, 2 ancillas, Heron, fractional": {
        "color": "tab:red",
        "marker": "v",
        "linestyle": "--",
    },
}

def plot_2q_depth_and_counts(
    sizes, results, xlabel="sites L", depth_ticks=(5, 10, 20, 50, 100)
):
    """Two-qubit depth (left, log scale) and two-qubit gate count (right) against the system size: one line per entry of `results`.

    `results` maps a label to one `Stats` per size (as returned by `route`).
    The shaded band reaches from the best to the worst of the seeds.
    Pure plotting: nothing is computed here. Returns the figure, so that a cell can end with the call."""
    figure, axes = plt.subplots(1, 2, figsize=(11, 4), layout="constrained")
    for label, series in results.items():
        style = STYLES[label]
        axes[0].plot(sizes, [s.depth for s in series], label=label, **style)
        axes[0].fill_between(
            sizes,
            [s.depth for s in series],
            [s.high for s in series],
            alpha=0.15,
            lw=0,
            color=style["color"],
        )
        axes[1].plot(sizes, [s.count for s in series], label=label, **style)
    axes[0].set(
        xlabel=xlabel,
        ylabel=f"two-qubit depth of the best routed circuit\n(band up to the worst of {NUM_SEEDS} seeds)",
        title="depth",
        yscale="log",
    )
    axes[0].yaxis.set_major_formatter(plt.ScalarFormatter())
    axes[0].set_yticks(list(depth_ticks))
    axes[1].set(
        xlabel=xlabel,
        ylabel="two-qubit gates of that circuit",
        title="two-qubit gates",
    )
    for axis in axes:
        axis.grid(axis="y")
    axes[1].legend(fontsize=8)
    plt.close(
        figure
    )  # drawn once, by being the cell's value; otherwise the inline backend shows it a second time
    return figure

def offset_in_L(offset, num_sites):
    """Show an offset in terms of the number of sites L: constants stay numbers, offsets that grow with L read `L + c`.

    This is a display rule, not a fit. It assumes every constant offset in the tables below is smaller than L (the largest
    is 3, from `ZXXZ`) and that no offset reaches 2L, which would have to be written `L + L`. The asserts in `pattern_table`
    guard that, so a size where the rule would mislabel an offset fails loudly instead.
    """
    if offset < num_sites:
        return str(offset)
    return "L" if offset == num_sites else f"L + {offset - num_sites}"


def pattern_table(operator, encoding, num_qubits, num_sites):
    """Per flow set: the Pauli patterns its terms take after encoding.

    A pattern is a Pauli string together with the offsets of the qubits it touches, measured from the first one:
    `ZZ @ (0, 2)` is a Z on some qubit q and a Z on qubit q + 2. This says which qubits of the *encoding register* a term
    couples. It says nothing about distance on a chip; that depends on where the layout puts those qubits. Offsets that
    grow with the number of sites `num_sites` (= L) are shown as `L + c`, so the scaling reads off directly.
    """
    assert (
        num_sites >= 4
    ), "constant offsets (up to 3) would be mistaken for offsets that grow with L"
    for index, group in enumerate(operator.split_out_groups()):
        encoded = [
            (label, idx)
            for label, idx, _ in encoding(group, num_qubits).to_sparse_list()
            if idx
        ]
        patterns = sorted(
            {
                (label, tuple(i - idx[0] for i in idx))
                for label, idx in encoded
            }
        )
        assert all(
            max(offsets) < 2 * num_sites for _, offsets in patterns
        ), "an offset reaches 2L"
        shown = ", ".join(
            f"{label} @ ({', '.join(offset_in_L(o, num_sites) for o in offsets)})"
            for label, offsets in patterns
        )
        print(
            f"  flow set {index}: {len(encoded):2d} terms, Pauli patterns: {shown}"
        )

def draw_flow_set(num_sites=L_SMALL):
    """The vertex and transfer operators of a spinless chain: the sites on a line, the east and west arrows as arcs below and above it.

    `mpl_draw` puts every edge label on the straight line between its two nodes, so the labels of the east and the west arrow
    of a bond would collide. We draw the arrows with `mpl_draw` and write the labels ourselves."""
    graph = rx.PyDiGraph()
    graph.add_nodes_from(range(num_sites))
    graph.add_edges_from([(i, i + 1, None) for i in range(num_sites - 1)])
    graph.add_edges_from([(i + 1, i, None) for i in range(num_sites - 1)])
    figure, axis = plt.subplots(
        figsize=(1.8 * num_sites, 2.2), layout="constrained"
    )
    for i in range(
        num_sites
    ):  # the grid of `draw_orderings`: the position of every site on the line
        axis.axvline(i, color="0.9", lw=1, zorder=0)
    axis.axhline(0, color="0.9", lw=1, zorder=0)
    mpl_draw(
        graph,
        pos={i: (i, 0) for i in range(num_sites)},
        ax=axis,
        labels=lambda v: f"$V_{v}$",
        with_labels=True,
        node_color=SPINLESS,
        node_size=600,
        edge_color=[HOP_EAST] * (num_sites - 1)
        + [HOP_WEST] * (num_sites - 1),
        width=2,
        connectionstyle="arc3,rad=0.3",
    )
    for i in range(num_sites - 1):
        axis.text(
            i + 0.5,
            -0.12,
            f"$T_{{{i},{i + 1}}}$",
            ha="center",
            va="center",
            fontsize=9,
            bbox=dict(
                boxstyle="round,pad=0.1", fc="white", ec="none", alpha=0.85
            ),
            zorder=5,
        )
        axis.text(
            i + 0.5,
            0.12,
            f"$T_{{{i + 1},{i}}}$",
            ha="center",
            va="center",
            fontsize=9,
            bbox=dict(
                boxstyle="round,pad=0.1", fc="white", ec="none", alpha=0.85
            ),
            zorder=5,
        )
    axis.set(xlim=(-0.5, num_sites - 0.5), ylim=(-0.4, 0.4))
    axis.axis("off")
    plt.close(
        figure
    )  # drawn once, by being the cell's value; otherwise the inline backend shows it a second time
    return figure


def draw_orderings(num_sites=L_SMALL):
    """Both mode orderings as two parallel chains (spin up on top, spin down below), each node at x = its mode index.

    The grid marks the position of every mode on the line. Hops are the arrows along a chain, the on-site interaction is the
    dotted on-site link between the two spins of a site. The ordering shows in the on-site links."""
    figure, axes = plt.subplots(2, 1, figsize=(9, 5.2), layout="constrained")
    for axis, order in zip(axes, ("blocked", "interleaved")):
        up, down = mode_maps(num_sites, order)
        pos = {
            mode(j): (mode(j), height)
            for mode, height in ((up, 0.5), (down, 0))
            for j in range(num_sites)
        }
        hops = [
            pair
            for mode in (up, down)
            for j in range(num_sites - 1)
            for pair in ((mode(j), mode(j + 1)), (mode(j + 1), mode(j)))
        ]
        graph = rx.PyDiGraph()
        graph.add_nodes_from(range(2 * num_sites))
        graph.add_edges_from([(a, b, None) for a, b in hops])
        for m in range(2 * num_sites):
            axis.axvline(m, color="0.9", lw=1, zorder=0)
        for height in (0, 0.5):
            axis.axhline(height, color="0.9", lw=1, zorder=0)
        mpl_draw(
            graph,
            pos=pos,
            ax=axis,
            labels=lambda m: f"$V_{{{m}}}$",
            with_labels=True,
            node_color=[
                SPIN_UP if pos[m][1] else SPIN_DOWN
                for m in range(2 * num_sites)
            ],
            node_size=600,
            edge_color=[HOP_EAST if a < b else HOP_WEST for a, b in hops],
            width=2,
            connectionstyle="arc3,rad=0.25",
        )
        for j in range(num_sites):
            axis.plot(
                *zip(pos[up(j)], pos[down(j)]),
                color=SPINLESS,
                lw=2,
                ls=":",
                zorder=1,
            )
        for (
            a,
            b,
        ) in hops:  # the transfer operator T_ab of every arrow, written next to the middle of its arc
            (xa, ya), (xb, _) = pos[a], pos[b]
            axis.text(
                (xa + xb) / 2,
                ya + (-0.12 if a < b else 0.12) * (1 + 0.1 * abs(b - a)),
                f"$T_{{{a},{b}}}$",
                ha="center",
                va="center",
                fontsize=9,
                bbox=dict(
                    boxstyle="round,pad=0.1",
                    fc="white",
                    ec="none",
                    alpha=0.85,
                ),
                zorder=5,
            )
        axis.text(-0.6, 0.5, "spin ↑", ha="right", va="center")
        axis.text(-0.6, 0, "spin ↓", ha="right", va="center")
        axis.set(
            title=f"{order} ordering",
            xlim=(-1.6, 2 * num_sites - 0.4),
            ylim=(-0.45, 0.95),
        )
        axis.axis("off")
    plt.close(
        figure
    )  # drawn once, by being the cell's value; otherwise the inline backend shows it a second time
    return figure


def ladder_strips(num_sites, orientation):
    """Strips needed for a chain of L + 1 qubits: it runs along 10 columns (horizontal) or 12 rows (vertical)."""
    return -(
        -(num_sites + 1) // (WIDTH if orientation == "horizontal" else ROWS)
    )


def ladder_orientation(num_sites):
    """The orientation with fewer strips (the vertical one never needs more), horizontal on a tie, among those that fit."""
    fitting = [
        o
        for o in ("horizontal", "vertical")
        if 2 * ladder_strips(num_sites, o)
        <= (ROWS if o == "horizontal" else WIDTH)
    ]
    if not fitting:
        raise ValueError(
            f"a ladder of {num_sites} sites does not fit on the {ROWS} x {WIDTH} lattice"
        )
    return min(
        fitting, key=lambda o: (ladder_strips(num_sites, o), o == "vertical")
    )


def folded_ladder(num_sites, orientation=None):
    """Virtual qubit -> physical qubit for the two-chain encoding.

    A strip is two adjacent grid lines, one per chain, so rung q joins the two lines at the same position along the
    strip. Strips are entered from alternating ends (a boustrophedon). Every strip boundary stretches one up link and
    one down link to distance 2, so the fewer the strips the fewer the stretched links.

    Horizontal strips run along the 10 columns and stack down the 12 rows; vertical strips run along the 12 rows and stack
    across the 10 columns. With `orientation=None` we take `ladder_orientation(num_sites)`."""
    orientation = orientation or ladder_orientation(num_sites)
    along_length = WIDTH if orientation == "horizontal" else ROWS
    layout = [None] * (2 * num_sites + 2)
    for q in range(num_sites + 1):
        strip, c = divmod(q, along_length)
        along = c if strip % 2 == 0 else along_length - 1 - c
        for chain, cross in ((0, 2 * strip), (num_sites + 1, 2 * strip + 1)):
            layout[chain + q] = (
                cross * WIDTH + along
                if orientation == "horizontal"
                else along * WIDTH + cross
            )
    assert (
        None not in layout
        and len(set(layout)) == len(layout)
        and max(layout) < nighthawk.size()
    )
    return layout


def draw_on_chip(
    axis, coupling_map, positions, layout, links, colours, labels, title
):
    """One layout drawn on the whole coupling map of a chip (unused qubits and couplers in gray).

    `positions[p]` is the (x, y) place of physical qubit `p` in the drawing and `layout[v]` the physical qubit of virtual
    qubit `v`. `links` are pairs of virtual qubits that the circuit couples along a chain or a lattice bond: a solid line if
    the pair is a coupler, a dotted red one if it is not. Each circle has the colour and the label given for its virtual qubit."""
    edges = {tuple(sorted(e)) for e in coupling_map.get_edges()}
    for a, b in edges:
        axis.plot(
            *zip(positions[a], positions[b]), color="0.85", lw=1, zorder=1
        )
    axis.scatter(
        *zip(*positions.values()),
        s=60,
        color="0.9",
        edgecolor="0.7",
        linewidth=0.6,
        zorder=2,
    )
    for q, r in links:
        stretched = tuple(sorted((layout[q], layout[r]))) not in edges
        start, end = positions[layout[q]], positions[layout[r]]
        style = dict(
            color="tab:red" if stretched else "0.4",
            lw=2 if stretched else 1.5,
            ls=":" if stretched else "-",
            zorder=3,
        )
        # a link that runs over the circle of another qubit of the layout is drawn as an arc, so that it is not mistaken for
        # a link that is missing between that qubit and its neighbors
        direction = np.subtract(end, start)
        covered = any(
            np.isclose(direction[0] * offset[1] - direction[1] * offset[0], 0)
            and 0 < np.dot(offset, direction) < np.dot(direction, direction)
            for offset in (
                np.subtract(positions[p], start)
                for p in layout
                if p not in (layout[q], layout[r])
            )
        )
        if covered:
            axis.add_patch(
                FancyArrowPatch(
                    start,
                    end,
                    connectionstyle="arc3,rad=-0.5",
                    arrowstyle="-",
                    **{**style, "linestyle": style.pop("ls")},
                )
            )
        else:
            axis.plot(*zip(start, end), **style)
    for virtual, physical in enumerate(layout):
        x, y = positions[physical]
        axis.scatter(x, y, s=330, color=colours[virtual], zorder=4)
        light = (
            sum(to_rgb(colours[virtual])) / 3 > 0.62
        )  # dark text on light circles
        axis.annotate(
            labels[virtual],
            (x, y),
            ha="center",
            va="center",
            fontsize=7,
            color="black" if light else "white",
            zorder=5,
        )
    xs, ys = zip(*positions.values())
    axis.set_xlim(
        min(xs) - 0.7, max(xs) + 0.7
    )  # room for the circles at the edge of the chip
    axis.set_ylim(min(ys) - 0.7, max(ys) + 0.7)
    axis.set_aspect("equal")
    axis.set_title(title, fontsize=9)
    axis.axis("off")


def chip_figure(num_panels):
    return plt.subplots(
        1,
        num_panels,
        figsize=(num_panels * (0.5 * WIDTH + 0.6), 0.5 * ROWS + 1.2),
        layout="constrained",
    )


def draw_fold(sizes):
    """The ladder layout for each size in `sizes`, side by side on the whole Nighthawk coupling map.

    Each circle carries the number of its virtual qubit, as in the drawings of section 2. The two ancillas are green and
    labeled `a↑` and `a↓`. A dotted red link is a chain link that is not a coupler."""
    figure, axes_fold = chip_figure(len(sizes))
    for axis, num_sites in zip(np.atleast_1d(axes_fold), sizes):
        up, down = (
            range(0, num_sites + 1),
            range(num_sites + 1, 2 * num_sites + 2),
        )
        ancillas = (num_sites, 2 * num_sites + 1)
        colours = [
            ANCILLA
            if v in ancillas
            else (SPIN_UP if v < num_sites else SPIN_DOWN)
            for v in range(2 * num_sites + 2)
        ]
        labels = [
            ("a↑" if v == num_sites else "a↓") if v in ancillas else str(v)
            for v in range(2 * num_sites + 2)
        ]
        links = [
            pair
            for chain in (up, down)
            for pair in zip(chain[:-1], chain[1:])
        ]
        orientation = ladder_orientation(num_sites)
        strips = ladder_strips(num_sites, orientation)
        title = (
            f"L = {num_sites}: {strips} {orientation} strip{'s' if strips > 1 else ''}\n"
            "up (dark blue), down (light blue), ancillas (green); dotted red = stretched link"
        )
        draw_on_chip(
            axis,
            nighthawk,
            nighthawk_positions,
            folded_ladder(num_sites),
            links,
            colours,
            labels,
            title,
        )
    plt.close(figure)  # shown once, as the cell's value
    return figure


def draw_interleaved(sizes):
    """The interleaved one-ancilla layouts for each size in `sizes` (one row each), Nighthawk and Heron side by side.

    Virtual qubit $2j$ is spin up and $2j + 1$ spin down on site $j$ (dark and light blue), the last one is the ancilla
    (green). The links are the chain links $q - (q + 1)$ of the encoded circuit."""
    # the layouts of the routing studies (section 4): a snake through the rows, a path found by search
    devices = (
        ("Nighthawk", nighthawk, nighthawk_positions, nighthawk_snake),
        (
            "Heron",
            heron,
            heron_positions,
            lambda num_qubits: find_path(heron, num_qubits),
        ),
    )
    figure, axes = plt.subplots(
        len(sizes),
        2,
        figsize=(0.4 * (WIDTH + 17) + 0.6, 0.4 * 16 * len(sizes) + 1.0),
        width_ratios=[WIDTH, 17],
        layout="constrained",
        squeeze=False,
    )
    for row, num_sites in zip(axes, sizes):
        num_qubits = 2 * num_sites + 1
        colours = [
            ANCILLA
            if v == 2 * num_sites
            else (SPIN_UP if v % 2 == 0 else SPIN_DOWN)
            for v in range(num_qubits)
        ]
        labels = [
            "a" if v == 2 * num_sites else str(v) for v in range(num_qubits)
        ]
        links = list(zip(range(num_qubits - 1), range(1, num_qubits)))
        for axis, (name, coupling_map, positions, layout_for) in zip(
            row, devices
        ):
            if name == "Nighthawk" and num_qubits > 120:
                print(
                    f"Cannot lay out {num_qubits} qubits on a 120-qubit Nighthawk chip!"
                )
                continue
            if name == "Heron" and num_qubits > 127:
                print(
                    f"Cannot lay out {num_qubits} qubits on a 1D chain on a Heron chip!"
                )
                continue
            layout = layout_for(num_qubits)
            title = f"{name}, L = {num_sites}: {num_qubits} qubits\nspin ↑ (dark blue), spin ↓ (light blue), ancilla (green)"
            draw_on_chip(
                axis,
                coupling_map,
                positions,
                layout,
                links,
                colours,
                labels,
                title,
            )
    plt.close(figure)  # shown once, as the cell's value
    return figure


def stripe_layout(rows, cols):
    """Virtual qubit -> Nighthawk qubit. Lattice column c is the chip row pair (2c, 2c + 1) and lattice row r is chip column r,
    so a lattice row is a stripe that alternates site, ancilla, site, ancilla, ... down one chip column."""
    if cols > ROWS // 2 or rows > WIDTH:
        raise ValueError(
            f"a {rows} x {cols} lattice does not fit the stripe layout on the {ROWS} x {WIDTH} chip"
        )
    layout = [None] * (2 * rows * cols)
    for r in range(rows):
        for c in range(cols):
            layout[r * cols + c] = (2 * c) * WIDTH + r  # the site
            layout[rows * cols + r * cols + c] = (
                2 * c + 1
            ) * WIDTH + r  # its ancilla
    return layout


def draw_stripes(shapes):
    """The stripe layout for each (rows, cols) in `shapes`, side by side on the whole Nighthawk coupling map.

    Circles carry the number of their virtual qubit: sites are $0 .. N-1$ (blue, row by row), their ancillas
    $N .. 2N-1$ (green). A solid line is a lattice bond between two sites on a coupler, a dotted red one is a bond without."""
    figure, axes_stripes = chip_figure(len(shapes))
    for axis, (rows, cols) in zip(np.atleast_1d(axes_stripes), shapes):
        num_sites = rows * cols
        links = [
            (j, j + 1) for j in range(num_sites) if j % cols < cols - 1
        ] + [(j, j + cols) for j in range(num_sites - cols)]
        colours = [SPINLESS] * num_sites + [ANCILLA] * num_sites
        labels = [str(v) for v in range(2 * num_sites)]
        title = f"{rows} x {cols} lattice, {2 * num_sites} qubits\nsites (blue), ancillas (green); dotted red = bond without a coupler"
        draw_on_chip(
            axis,
            nighthawk,
            nighthawk_positions,
            stripe_layout(rows, cols),
            links,
            colours,
            labels,
            title,
        )
    plt.close(figure)  # shown once, as the cell's value
    return figure

draw_flow_set()

def flow_set_groups(operator, stride=1):
    """0: arrows pointing 'east', 1: arrows pointing 'west', 2: the diagonal terms."""
    groups = []
    for terms, _ in operator.iter_terms():
        match terms:
            case [(left, right)] if right == left + stride:
                groups.append(0)
            case [(left, right)] if right == left - stride:
                groups.append(1)
            case _:
                groups.append(2)
    return groups


def spinless_chain(num_sites, tunneling, interaction):
    data = defaultdict(complex)
    for j in range(num_sites - 1):
        k = j + 1
        data[((j, k),)] += tunneling  # hopping: t (T_jk + T_kj)
        data[((k, j),)] += tunneling
        data[()] += (
            interaction / 4
        )  # interaction: U n_j n_k, with n_j = (1 - V_j) / 2
        data[((j, j),)] -= interaction / 4
        data[((k, k),)] -= interaction / 4
        data[((j, j), (k, k))] += interaction / 4
    operator = TransferVertexOperator.from_dict(data)
    operator.groups = flow_set_groups(operator)
    return operator


spinless = spinless_chain(L_SMALL, 1.0, 2.0)
print(f"spinless chain, {L_SMALL} sites")
for group, op in enumerate(spinless.split_out_groups()):
    print(f"{group=}\n{format(op)}")

def parity_delocalized(action, num_qubits):
    """Image of ONE vertex/transfer operator under the one-ancilla parity-delocalized encoding.

    The guide needs the first three cases only. The last two (modes TWO apart) appear once spins are interleaved.
    """
    match action:
        case (j, k) if j == k:  # V_j = Z_j Z_{j+1}
            return SparseObservable.from_sparse_list(
                [("ZZ", [j, j + 1], 1.0)], num_qubits=num_qubits
            )
        case (j, k) if k == j + 1:  # weight 1
            # east flow-set in spinless or blocked spinful order
            return SparseObservable.from_sparse_list(
                [("X", [k], -0.5)], num_qubits=num_qubits
            )
        case (j, k) if k == j - 1:  # weight 3
            # west flow-set in spinless or blocked spinful order
            return SparseObservable.from_sparse_list(
                [("ZXZ", [k, k + 1, k + 2], 0.5)], num_qubits=num_qubits
            )
        case (j, k) if k == j + 2:  # weight 2
            # east flow-set in interleaved spinful order
            return SparseObservable.from_sparse_list(
                [("YY", [j + 1, j + 2], 0.5)], num_qubits=num_qubits
            )
        case (j, k) if k == j - 2:  # weight 4
            # west flow-set in interleaved spinful order
            return SparseObservable.from_sparse_list(
                [("ZXXZ", [k, k + 1, k + 2, k + 3], 0.5)],
                num_qubits=num_qubits,
            )
        case _:
            raise ValueError(f"no image for transfer operator {action}")


def encode(operator, num_qubits, action=parity_delocalized):
    """Map a whole operator by composing the images of its generators."""
    return map_transfer_vertex_generators(
        operator,
        lambda generator: action(generator, num_qubits),
        identity=lambda: SparseObservable.identity(num_qubits),
        compose=SparseObservable.compose,
    ).simplify()


print(f"flow-set encoded spinless chain -> {L_SMALL + 1} qubits")
pattern_table(spinless, encode, L_SMALL + 1, L_SMALL)