Are there qiskit or stim examples of how to decode a subsystem code

Question

I can manage decoding a stabilizer code in qiskit and stim. I'd like to try a subsystem code next. For a stabilizer code, all stabilizers commute so simultaneous measurement of these produces predictable results. For subsystem codes (such as bacon-shor, bacon-bravyi, ...) you measure (gauge) operators that do not necessarily commute. I'd like to see how this is handled in qiskit or stim; also how the gauge operator measurements are assembled into syndromes of the stabilizer subgroup.

Answer 1

The main concept you need to know, to encode a subsystem code into Stim, is that Stim circuits have no concept of a "code". But Stim circuits understand producing detection events by xoring measurements together. And you can explain a subsystem code in that language.

The paper "Less Bacon More Threshold" includes working code on Zenodo. (There's also a version on github.) The code's not really intended to be pedagogical, but it works.

In the case of the Bacon-Shor code, what this looks like is that when you declare a DETECTOR you need to give it an entire column (or row) of measurements from two adjacent rounds. That encodes the comparison of the two stabilizers, which are formed in each round by multiplying together those measurements. In the code one implementation of that looks like this:

def _det(
    *,
    builder: gen.Builder,
    group: List[gen.Tile],
    layers: List[Any],
    include_data: bool = False,
):
    builder.detector(
        [
            gen.AtLayer(q, layer)
            for tile in group
            for layer in layers
            for q in (tile.used_set if include_data else [tile.measurement_qubit])
        ],
        pos=min([tile.measurement_qubit for tile in group], key=gen.complex_key),
    )


def _do_det_groups(
    *,
    builder: gen.Builder,
    patch: gen.Patch,
    basis: str,
    layers: List[Any],
    include_data: bool = False
):
    groups = sinter.group_by(
        [tile for tile in patch.tiles if tile.basis == basis],
        key=lambda tile:
            tile.measurement_qubit.real
            if basis == 'X'
            else tile.measurement_qubit.imag,
    )
    for _, group in sorted(groups.items()):
        _det(builder=builder, group=group, layers=layers, include_data=include_data)

That code is using various utilities internal to that code, but you can at least see that it's grouping into rows/columns and then merging measurements along those groups.

In the actual stim circuit it looks like this:

QUBIT_COORDS(0, 0) 0
QUBIT_COORDS(0, 1) 1
QUBIT_COORDS(0, 2) 2
QUBIT_COORDS(0, 3) 3
QUBIT_COORDS(1, 0) 4
QUBIT_COORDS(1, 1) 5
QUBIT_COORDS(1, 2) 6
QUBIT_COORDS(1, 3) 7
QUBIT_COORDS(2, 0) 8
QUBIT_COORDS(2, 1) 9
QUBIT_COORDS(2, 2) 10
QUBIT_COORDS(2, 3) 11
QUBIT_COORDS(3, 0) 12
QUBIT_COORDS(3, 1) 13
QUBIT_COORDS(3, 2) 14
QUBIT_COORDS(3, 3) 15
R 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
TICK
MPP X0*X4 X1*X5 X2*X6 X3*X7 X8*X12 X9*X13 X10*X14 X11*X15
TICK
MPP X4*X8 X5*X9 X6*X10 X7*X11
TICK
MPP Z0*Z1 Z2*Z3 Z4*Z5 Z6*Z7 Z8*Z9 Z10*Z11 Z12*Z13 Z14*Z15
TICK
MPP Z1*Z2 Z5*Z6 Z9*Z10 Z13*Z14
DETECTOR(0, 0.5, 0) rec[-12] rec[-10] rec[-8] rec[-6]
DETECTOR(0, 1.5, 0) rec[-4] rec[-3] rec[-2] rec[-1]
DETECTOR(0, 2.5, 0) rec[-11] rec[-9] rec[-7] rec[-5]
SHIFT_COORDS(0, 0, 1)
TICK
REPEAT 98 {
    MPP X0*X4 X1*X5 X2*X6 X3*X7 X8*X12 X9*X13 X10*X14 X11*X15
    TICK
    MPP X4*X8 X5*X9 X6*X10 X7*X11
    TICK
    MPP Z0*Z1 Z2*Z3 Z4*Z5 Z6*Z7 Z8*Z9 Z10*Z11 Z12*Z13 Z14*Z15
    TICK
    MPP Z1*Z2 Z5*Z6 Z9*Z10 Z13*Z14
    DETECTOR(0.5, 0, 0) rec[-48] rec[-47] rec[-46] rec[-45] rec[-24] rec[-23] rec[-22] rec[-21]
    DETECTOR(1.5, 0, 0) rec[-40] rec[-39] rec[-38] rec[-37] rec[-16] rec[-15] rec[-14] rec[-13]
    DETECTOR(2.5, 0, 0) rec[-44] rec[-43] rec[-42] rec[-41] rec[-20] rec[-19] rec[-18] rec[-17]
    DETECTOR(0, 0.5, 0) rec[-36] rec[-34] rec[-32] rec[-30] rec[-12] rec[-10] rec[-8] rec[-6]
    DETECTOR(0, 1.5, 0) rec[-28] rec[-27] rec[-26] rec[-25] rec[-4] rec[-3] rec[-2] rec[-1]
    DETECTOR(0, 2.5, 0) rec[-35] rec[-33] rec[-31] rec[-29] rec[-11] rec[-9] rec[-7] rec[-5]
    SHIFT_COORDS(0, 0, 1)
    TICK
}
MPP X0*X4 X1*X5 X2*X6 X3*X7 X8*X12 X9*X13 X10*X14 X11*X15
TICK
MPP X4*X8 X5*X9 X6*X10 X7*X11
TICK
MPP Z0*Z1 Z2*Z3 Z4*Z5 Z6*Z7 Z8*Z9 Z10*Z11 Z12*Z13 Z14*Z15
TICK
MPP Z1*Z2 Z5*Z6 Z9*Z10 Z13*Z14
TICK
M 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
DETECTOR(0, 0.5, 0) rec[-28] rec[-16] rec[-15]
DETECTOR(0, 1.5, 0) rec[-20] rec[-15] rec[-14]
DETECTOR(0, 2.5, 0) rec[-27] rec[-14] rec[-13]
DETECTOR(1, 0.5, 0) rec[-26] rec[-12] rec[-11]
DETECTOR(1, 1.5, 0) rec[-19] rec[-11] rec[-10]
DETECTOR(1, 2.5, 0) rec[-25] rec[-10] rec[-9]
DETECTOR(2, 0.5, 0) rec[-24] rec[-8] rec[-7]
DETECTOR(2, 1.5, 0) rec[-18] rec[-7] rec[-6]
DETECTOR(2, 2.5, 0) rec[-23] rec[-6] rec[-5]
DETECTOR(3, 0.5, 0) rec[-22] rec[-4] rec[-3]
DETECTOR(3, 1.5, 0) rec[-17] rec[-3] rec[-2]
DETECTOR(3, 2.5, 0) rec[-21] rec[-2] rec[-1]
DETECTOR(0.5, 0, 0) rec[-64] rec[-63] rec[-62] rec[-61] rec[-40] rec[-39] rec[-38] rec[-37]
DETECTOR(1.5, 0, 0) rec[-56] rec[-55] rec[-54] rec[-53] rec[-32] rec[-31] rec[-30] rec[-29]
DETECTOR(2.5, 0, 0) rec[-60] rec[-59] rec[-58] rec[-57] rec[-36] rec[-35] rec[-34] rec[-33]
DETECTOR(0, 0.5, 0) rec[-52] rec[-50] rec[-48] rec[-46] rec[-28] rec[-26] rec[-24] rec[-22]
DETECTOR(0, 1.5, 0) rec[-44] rec[-43] rec[-42] rec[-41] rec[-20] rec[-19] rec[-18] rec[-17]
DETECTOR(0, 2.5, 0) rec[-51] rec[-49] rec[-47] rec[-45] rec[-27] rec[-25] rec[-23] rec[-21]
OBSERVABLE_INCLUDE(0) rec[-16] rec[-12] rec[-8] rec[-4]
SHIFT_COORDS(0, 0, 1)
TICK

The main distinction, compared to the stim circuit for a stabilizer code, is that the DETECTOR lines include a lot more targets.