Fix HERALDED_PAULI_CHANNEL_1 targeting the wrong qubits

src/stim/simulators/frame_simulator.inl

@@ -768,16 +768,17 @@ void FrameSimulator<W>::do_HERALDED_PAULI_CHANNEL_1(const CircuitInstruction &in
     RareErrorIterator::for_samples(t, nt * batch_size, rng, [&](size_t s) {
         auto shot = s % batch_size;
         auto target = s / batch_size;
+        auto qubit = inst.targets[target].qubit_value();
         m_record.storage[m_record.stored + target][shot] = 1;
 
         double p = dist(rng) * t;
         if (p < hx) {
-            x_table[target][shot] ^= 1;
+            x_table[qubit][shot] ^= 1;
         } else if (p < hx + hz) {
-            z_table[target][shot] ^= 1;
+            z_table[qubit][shot] ^= 1;
         } else if (p < hx + hz + hy) {
-            x_table[target][shot] ^= 1;
-            z_table[target][shot] ^= 1;
+            x_table[qubit][shot] ^= 1;
+            z_table[qubit][shot] ^= 1;
         }
     });
 

src/stim/simulators/frame_simulator.test.cc

@@ -1615,3 +1615,69 @@ TEST_EACH_WORD_SIZE_W(FrameSimulator, heralded_pauli_channel_1_statistics, {
     EXPECT_NEAR(bins[6] / (double)n, 0.25, 0.04);
     EXPECT_NEAR(bins[7] / (double)n, 0.15, 0.04);
 })
+
+TEST_EACH_WORD_SIZE_W(FrameSimulator, heralded_erase_statistics_offset_by_2, {
+    auto circuit = Circuit(R"CIRCUIT(
+        MXX 2 3
+        MZZ 2 3
+        HERALDED_ERASE(0.1) 2
+        MXX 2 3
+        MZZ 2 3
+        DETECTOR rec[-1] rec[-4]
+        DETECTOR rec[-2] rec[-5]
+        DETECTOR rec[-3]
+    )CIRCUIT");
+    size_t n;
+    std::array<size_t, 8> bins{};
+    FrameSimulator<W> sim(
+        circuit.compute_stats(), FrameSimulatorMode::STORE_DETECTIONS_TO_MEMORY, 1024, INDEPENDENT_TEST_RNG());
+    for (n = 0; n < 1024 * 256; n += 1024) {
+        sim.reset_all();
+        sim.do_circuit(circuit);
+        auto sample = sim.det_record.storage.transposed();
+        for (size_t k = 0; k < 1024; k++) {
+            bins[sample[k].u8[0]]++;
+        }
+    }
+    EXPECT_NEAR(bins[0] / (double)n, 0.9, 0.05);
+    EXPECT_EQ(bins[1], 0);
+    EXPECT_EQ(bins[2], 0);
+    EXPECT_EQ(bins[3], 0);
+    EXPECT_NEAR(bins[4] / (double)n, 0.025, 0.02);
+    EXPECT_NEAR(bins[5] / (double)n, 0.025, 0.02);
+    EXPECT_NEAR(bins[6] / (double)n, 0.025, 0.02);
+    EXPECT_NEAR(bins[7] / (double)n, 0.025, 0.02);
+})
+
+TEST_EACH_WORD_SIZE_W(FrameSimulator, heralded_pauli_channel_1_statistics_offset_by_2, {
+    auto circuit = Circuit(R"CIRCUIT(
+        MXX 2 3
+        MZZ 2 3
+        HERALDED_PAULI_CHANNEL_1(0.05, 0.10, 0.15, 0.25) 2
+        MXX 2 3
+        MZZ 2 3
+        DETECTOR rec[-1] rec[-4]
+        DETECTOR rec[-2] rec[-5]
+        DETECTOR rec[-3]
+    )CIRCUIT");
+    size_t n;
+    std::array<size_t, 8> bins{};
+    FrameSimulator<W> sim(
+        circuit.compute_stats(), FrameSimulatorMode::STORE_DETECTIONS_TO_MEMORY, 1024, INDEPENDENT_TEST_RNG());
+    for (n = 0; n < 1024 * 256; n += 1024) {
+        sim.reset_all();
+        sim.do_circuit(circuit);
+        auto sample = sim.det_record.storage.transposed();
+        for (size_t k = 0; k < 1024; k++) {
+            bins[sample[k].u8[0]]++;
+        }
+    }
+    EXPECT_NEAR(bins[0] / (double)n, 0.45, 0.05);
+    EXPECT_EQ(bins[1], 0);
+    EXPECT_EQ(bins[2], 0);
+    EXPECT_EQ(bins[3], 0);
+    EXPECT_NEAR(bins[4] / (double)n, 0.05, 0.04);
+    EXPECT_NEAR(bins[5] / (double)n, 0.10, 0.04);
+    EXPECT_NEAR(bins[6] / (double)n, 0.25, 0.04);
+    EXPECT_NEAR(bins[7] / (double)n, 0.15, 0.04);
+})

src/stim/simulators/frame_simulator_pybind_test.py

@@ -150,29 +150,6 @@ def test_peek_pauli_flips():
     assert 250 < v4[0] < 450
 
 
-def test_surface_code():
-    circuit = stim.Circuit.generated(
-        "surface_code:rotated_memory_x",
-        distance=3,
-        rounds=5,
-        before_round_data_depolarization=1e-2,
-        before_measure_flip_probability=1e-2,
-        after_reset_flip_probability=1e-2,
-        after_clifford_depolarization=1e-2,
-    )
-
-    # Find the index of the final MR layer.
-    mr_layer = len(circuit) - 1
-    while circuit[mr_layer].name != 'MR':
-        mr_layer -= 1
-    circuit_before_mr = circuit[:mr_layer]
-
-    sim = stim.FlipSimulator(batch_size=256, disable_stabilizer_randomization=True)
-    sim.do(circuit_before_mr)
-    for b in sim.peek_pauli_flips():
-        print(list(b))
-
-
 def test_set_pauli_flip():
     sim = stim.FlipSimulator(
         batch_size=2,
@@ -238,3 +215,14 @@ def test_set_pauli_flip():
         stim.PauliString('____X'),
         stim.PauliString('XZ___'),
     ]
+
+
+def test_repro_heralded_pauli_channel_1_bug():
+    circuit = stim.Circuit("""
+        R 0 1
+        HERALDED_PAULI_CHANNEL_1(0.2, 0.2, 0, 0) 1
+        M 0
+    """)
+    result = circuit.compile_sampler().sample(1024)
+    assert np.sum(result[:, 0]) > 0
+    assert np.sum(result[:, 1]) == 0