Fix order of observable annotations not being reflected in sampler output (#127)

## Summary
- Fix `parse_stim_circuit` crash on empty `DETECTOR` /
`OBSERVABLE_INCLUDE` annotations; they now produce deterministic-zero
detector/observable bits, matching Stim.
- Fix sparse and out-of-order `OBSERVABLE_INCLUDE` indices: sampler
output now has one column per id in `range(num_observables)`, in sorted
order, with missing ids as deterministic-zero columns.

---------

Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>

Author: rafaelha

CHANGELOG.md

@@ -9,6 +9,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 - `MR`, `MRX`, and `MRY` no longer double-count their measurement flip probability as both a pre-measurement Pauli error and a measurement-result flip.
+- Out-of-order `OBSERVABLE_INCLUDE` indices now produce the correct sampler column order and output shape. Missing indices below the maximum mentioned id appear as deterministic-zero columns, and columns are emitted in sorted logical-index order.
+- Empty `DETECTOR` and `OBSERVABLE_INCLUDE` annotations (without targets) no longer crash the parser; they now produce zero detector/observable bits, matching Stim semantics.
+
 
 ### Added
 - `TPP` and `TPP_DAG` instructions — applies exp(-i pi/8 P) or exp(+i pi/8 P) (up to global phase) for a Pauli product P, i.e., phases the -1 eigenspace of P by exp(i pi/4) or exp(-i pi/4).

src/tsim/core/graph.py

@@ -231,7 +231,7 @@ def build_sampling_graph(
             g.remove_vertex(vertices.pop())
 
         labels = [f"det[{i}]" for i in range(len(built.detectors))] + [
-            f"obs[{i}]" for i in built.observables_dict
+            f"obs[{i}]" for i in sorted(built.observables_dict)
         ]
         for label in labels:
             vs = annotation_to_vertex[label]

src/tsim/core/instructions.py

@@ -1114,12 +1114,26 @@ def ry(b: GraphRepresentation, qubit: int) -> None:
 # =============================================================================
 
 
-def detector(b: GraphRepresentation, rec: list[int], *args) -> None:
-    """Add detector annotation that XORs the given measurement record bits."""
-    row = min({b.graph.row(b.rec[r]) for r in rec}) - 0.5
+def _annotation_row(b: GraphRepresentation, rec: list[int]) -> float:
+    """Pick a fresh row for a detector/observable vertex.
+
+    Empty `rec` is valid Stim and represents a deterministic-zero annotation;
+    fall back to a row above existing annotations rather than `min()` on an
+    empty set.
+    """
     d_rows = {b.graph.row(d) for d in b.detectors + b.observables}
+    if rec:
+        row: float = min(b.graph.row(b.rec[r]) for r in rec) - 0.5
+    else:
+        row = (max(d_rows) + 1) if d_rows else 0
     while row in d_rows:
         row += 1
+    return row
+
+
+def detector(b: GraphRepresentation, rec: list[int], *args) -> None:
+    """Add detector annotation that XORs the given measurement record bits."""
+    row = _annotation_row(b, rec)
     v0 = b.graph.add_vertex(
         VertexType.X, qubit=-1, row=row, phase=f"det[{len(b.detectors)}]"
     )
@@ -1133,10 +1147,7 @@ def observable_include(b: GraphRepresentation, rec: list[int], idx: int) -> None
     idx = int(idx)
 
     if idx not in b.observables_dict:
-        row = min({b.graph.row(b.rec[r]) for r in rec}) - 0.5
-        d_rows = {b.graph.row(d) for d in b.detectors + b.observables}
-        while row in d_rows:
-            row += 1
+        row = _annotation_row(b, rec)
         v0 = b.graph.add_vertex(VertexType.X, qubit=-1, row=row, phase=f"obs[{idx}]")
         b.observables_dict[idx] = v0
 

src/tsim/core/parse.py

@@ -261,4 +261,13 @@ def parse_stim_circuit(
                 gate_func(b, *chunk, *args)
 
     finalize_correlated_error(b)
+
+    # Materialize every observable id from 0..num_observables-1 so missing
+    # indices appear as deterministic-zero outputs and downstream iteration
+    # is in sorted index order, matching Stim semantics.
+    for i in range(stim_circuit.num_observables):
+        if i not in b.observables_dict:
+            observable_include(b, [], i)
+    b.observables_dict = {i: b.observables_dict[i] for i in sorted(b.observables_dict)}
+
     return b

src/tsim/utils/diagram.py

@@ -499,9 +499,10 @@ def render_pyzx_d3(stim_circ: stim.Circuit, kwargs: dict[str, Any]) -> GraphS:
         return g
 
     g = g.clone()
-    max_row = max(g.row(v) for v in built.last_vertex.values())
-    for q in built.last_vertex:
-        g.set_row(built.last_vertex[q], max_row)
+    if built.last_vertex:
+        max_row = max(g.row(v) for v in built.last_vertex.values())
+        for q in built.last_vertex:
+            g.set_row(built.last_vertex[q], max_row)
 
     for v in list(g.vertices()):
         phase_vars = g._phaseVars[v]

test/unit/core/test_parse.py

@@ -697,6 +697,38 @@ def test_spp_multiple_products(self):
         assert len(b.rec) == 0
 
 
+class TestParseEmptyAnnotations:
+    """Empty DETECTOR / OBSERVABLE_INCLUDE annotations are valid Stim and represent
+    deterministic-zero detector/observable bits. They must parse without crashing
+    and produce a single annotation vertex with no record edges."""
+
+    def test_empty_detector_alone(self):
+        b = parse_stim_circuit(stim.Circuit("DETECTOR"))
+        assert len(b.detectors) == 1
+        v = b.detectors[0]
+        assert b.graph.type(v) == VertexType.X
+        assert list(b.graph.neighbors(v)) == []
+
+    def test_empty_observable_alone(self):
+        b = parse_stim_circuit(stim.Circuit("OBSERVABLE_INCLUDE(0)"))
+        assert set(b.observables_dict) == {0}
+        v = b.observables_dict[0]
+        assert b.graph.type(v) == VertexType.X
+        assert list(b.graph.neighbors(v)) == []
+
+    def test_empty_detector_after_measurement(self):
+        b = parse_stim_circuit(stim.Circuit("M 0\nDETECTOR rec[-1]\nDETECTOR"))
+        assert len(b.detectors) == 2
+        # First detector has the measurement edge, second has no edges.
+        assert len(list(b.graph.neighbors(b.detectors[0]))) == 1
+        assert list(b.graph.neighbors(b.detectors[1])) == []
+
+    def test_empty_detector_with_args(self):
+        """DETECTOR(coords...) with empty rec must also parse."""
+        b = parse_stim_circuit(stim.Circuit("DETECTOR(1, 2)"))
+        assert len(b.detectors) == 1
+
+
 class TestParseMPPCancellation:
     """Tests for MPP with duplicate/anticommuting Pauli targets.
 
@@ -853,3 +885,30 @@ def test_tpp_anti_hermitian_raises(self):
         """TPP Z0*X0 = iY is anti-Hermitian and should raise."""
         with pytest.raises(ValueError, match="anti-Hermitian"):
             parse_stim_circuit(stim.Circuit("SPP[T] Z0*X0"))
+
+
+class TestParseSparseObservables:
+    """OBSERVABLE_INCLUDE indices are sparse and out-of-order in real circuits.
+    Stim defines num_observables = max_index + 1 and emits one column per id
+    in sorted order, with missing ids as deterministic-zero columns."""
+
+    def test_sparse_observable_index_pads_missing(self):
+        circuit = stim.Circuit("M 0\nOBSERVABLE_INCLUDE(2) rec[-1]")
+        b = parse_stim_circuit(circuit)
+        assert set(b.observables_dict) == {0, 1, 2}
+        # Missing ids 0 and 1 have no record edges (deterministic zero).
+        assert list(b.graph.neighbors(b.observables_dict[0])) == []
+        assert list(b.graph.neighbors(b.observables_dict[1])) == []
+        # Index 2 has the measurement edge.
+        assert len(list(b.graph.neighbors(b.observables_dict[2]))) == 1
+
+    def test_observables_dict_is_sorted_after_out_of_order(self):
+        circuit = stim.Circuit(
+            "M 0\nM 1\nOBSERVABLE_INCLUDE(2) rec[-2]\nOBSERVABLE_INCLUDE(0) rec[-1]"
+        )
+        b = parse_stim_circuit(circuit)
+        assert list(b.observables_dict) == [0, 1, 2]
+
+    def test_no_observables_remains_empty(self):
+        b = parse_stim_circuit(stim.Circuit("M 0"))
+        assert b.observables_dict == {}

test/unit/test_sampler.py

@@ -363,3 +363,63 @@ def test_reference_sample_defaults_unchanged():
         use_observable_reference_sample=False,
     )
     assert np.array_equal(d1, d2)
+
+
+def test_detector_sampler_empty_annotations():
+    """Empty DETECTOR / OBSERVABLE_INCLUDE produce deterministic-zero columns
+    that sit alongside non-empty ones, matching Stim semantics."""
+    c = Circuit("""
+        X 0
+        M 0 1
+        DETECTOR rec[-2]
+        DETECTOR
+        OBSERVABLE_INCLUDE(0) rec[-1]
+        OBSERVABLE_INCLUDE(1)
+        """)
+    assert c.num_detectors == 2
+    assert c.num_observables == 2
+
+    sampler = c.compile_detector_sampler(seed=0)
+    det, obs = sampler.sample(4, separate_observables=True)
+    assert det.shape == (4, 2)
+    assert obs.shape == (4, 2)
+    # First detector reads the X-induced 1; second detector and both observables
+    # are deterministically zero.
+    assert np.all(det[:, 0])
+    assert not np.any(det[:, 1])
+    assert not np.any(obs)
+
+
+def test_detector_sampler_sparse_observable_index():
+    """OBSERVABLE_INCLUDE(2) only must produce 3 observable columns, with the
+    measured value landing in column 2 and columns 0/1 deterministically zero."""
+    c = Circuit("X 0\nM 0\nOBSERVABLE_INCLUDE(2) rec[-1]")
+    assert c.num_observables == 3
+
+    sampler = c.compile_detector_sampler(seed=0)
+    samples = sampler.sample(4, append_observables=True)
+    assert samples.shape == (4, 3)
+    assert not np.any(samples[:, 0])
+    assert not np.any(samples[:, 1])
+    assert np.all(samples[:, 2])
+
+
+def test_detector_sampler_out_of_order_observable_indices():
+    """Out-of-order OBSERVABLE_INCLUDE must produce columns in sorted index order."""
+    # rec[-2] is qubit-0 (=1 after X), rec[-1] is qubit-1 (=0).
+    # Map: obs[2] = rec[-2] = 1, obs[0] = rec[-1] = 0, obs[1] is unmentioned = 0.
+    c = Circuit("""
+        X 0
+        M 0 1
+        OBSERVABLE_INCLUDE(2) rec[-2]
+        OBSERVABLE_INCLUDE(0) rec[-1]
+        """)
+    assert c.num_observables == 3
+
+    sampler = c.compile_detector_sampler(seed=0)
+    _, obs = sampler.sample(2, separate_observables=True)
+    assert obs.shape == (2, 3)
+    # Sorted-by-index order: [obs0, obs1, obs2] = [0, 0, 1]
+    assert not np.any(obs[:, 0])
+    assert not np.any(obs[:, 1])
+    assert np.all(obs[:, 2])