How to find a set of independent logical operators for a stabilizer code with Stim

As a simple example, say I am studying the stabilizer code with $$n=4$$ physical qubits and stabilizers $$Z_0Z_1Z_2Z_3$$ and $$X_0X_1X_2X_3$$ (I'll be using 0-based indexing throughout). Using Stim, I could get a representation of the code space by initializing $$n$$ Bell pairs and measuring the stabilizers on one side:

import stim

s = stim.TableauSimulator()
s.h(0, 1, 2, 3)
s.cx(0, 4, 1, 5, 2, 6, 3, 7)
s.measure_observable(stim.PauliString('ZZZZ'))
s.measure_observable(stim.PauliString('XXXX'))

Now s.canonical_stabilizers() gives:

[stim.PauliString("+X__X_XX_"),
stim.PauliString("+Z__Z_ZZ_"),
stim.PauliString("+_X_X_X_X"),
stim.PauliString("+_Z_Z_Z_Z"),
stim.PauliString("+__XX__XX"),
stim.PauliString("+__ZZ__ZZ"),
stim.PauliString("+____XXXX"),
stim.PauliString("+____ZZZZ")]

(The signs depend on random measurement results, and for now I don't care about them.)

I know that this stabilizer group has the following structure: Each stabilizer with support only on one "side" (either qubits 0, 1, 2, 3 or qubits 4, 5, 6, 7) corresponds to a code stabilizer (e.g. ZZZZ____, XXXX____, ____ZZZZ, ____XXXX), while each stabilizer with support on both sides and cannot be written as a combination of the previous kind of stabilizers corresponds to a logical operator. For example, the first 6 lines of s.canonical_stabilizers() suggest that $$X_0X_3$$, $$Z_0Z_3$$, $$X_1X_3$$, $$Z_1Z_3$$, $$X_2X_3$$, and $$Z_2Z_3$$ are logical operators of the code.

Now, this is all good, but these logical operators are not all independent up to stabilizers: For example, $$X_0X_3 \cdot X_1X_3 \cdot X_2X_3 = X_0X_1X_2X_3$$ is a stabilizer. If I could find a set of generators of the abovementioned 8-qubit stabilizer group that includes all of the stabilizers on both sides (i.e. ZZZZ____, XXXX____, ____ZZZZ, ____XXXX in this case), then the remaining stabilizers will give me a set of independent logical operators.

I know that I could make this happen with some sort of Gaussian elimination procedure, but an actually implementation of that would be cumbersome, error-prone, and maybe not very fast if using Python. Is there an idiomatic way in Stim of finding out a set of logical operators? (It would be nice if they satisfy the canonical commutation relations too, but currently I have no idea how to do this even not considering the commutation relations.)

stim.Tableau.from_stabilizers will solve for the observables as part of completing a tableau. It works by finding operations that turn the stabilizers into single-qubit terms. The observables are then created by looking at what undoing those operations turns the other qubits into.

import stim

def get_observables(
stabilizers: list[stim.PauliString],
) -> list[tuple[stim.PauliString, stim.PauliString]]:
completed_tableau = stim.Tableau.from_stabilizers(
stabilizers,
allow_redundant=True,
allow_underconstrained=True,
)

observables = []
for k in range(len(completed_tableau))[::-1]:
z = completed_tableau.z_output(k)
if z in stabilizers:
break
x = completed_tableau.x_output(k)
observables.append((x, z))

return observables

print("[[4,2,2]] iceberg code observables")
for obs in get_observables([
stim.PauliString("XXXX"),
stim.PauliString("ZZZZ"),
]):
print(obs)

# prints:
# [[4,2,2]] iceberg code observables
# (stim.PauliString("+_X_X"), stim.PauliString("+Z__Z"))
# (stim.PauliString("+_XX_"), stim.PauliString("+Z_Z_"))

print("[[5,1,3]] perfect code observables")
for obs in get_observables([
stim.PauliString("XZZX_"),
stim.PauliString("_XZZX"),
stim.PauliString("X_XZZ"),
stim.PauliString("ZX_XZ"),
]):
print(obs)
# prints:
# [[5,1,3]] perfect code observables
# (stim.PauliString("-Z_XX_"), stim.PauliString("-_ZXZ_"))