Improve STN/MAST simulators: locality-ordered projection, prefix-tree sampling, qec-grade defaults, full Python exposure and docs

docs/user-guide/stabilizer-tensor-networks.md

@@ -0,0 +1,145 @@
+# Stabilizer Tensor-Network Simulators
+
+The experimental `pecos_rslib_exp` package exposes three related tools:
+
+- `StabMps` executes Clifford and arbitrary-rotation circuits using a stabilizer tableau plus an MPS of coefficients.
+- `Mast` defers non-Clifford work through preallocated magic-state ancillas.
+- `StabMpsCompile` replays a circuit without an MPS and estimates which execution strategy is practical.
+
+All public bitstrings use qubit-index order: `bits[q]` is the bit for qubit `q`. Dense state-vector indices are little-endian, so a row maps to `sum(int(bits[q]) << q for q in range(len(bits)))`. Gate rotation angles are in radians.
+
+## `StabMps` quickstart
+
+Use `sample_bitstrings`, plural, for shot workloads. It shares each distinct measurement-prefix projection between shots; `sample_bitstring` clones and collapses the entire simulator once per shot.
+
+```python
+import math
+
+import pecos_rslib_exp as exp
+
+sim = exp.StabMps(2, seed=7, lazy_measure=True)
+sim.run_gate("H", {0})
+sim.run_gate("CX", {(0, 1)})
+sim.run_gate("RZ", {1}, angle=math.pi / 4)
+
+shots = sim.sample_bitstrings(32)
+assert len(shots) == 32
+assert all(len(bits) == 2 for bits in shots)
+assert all(bits in ([False, False], [True, True]) for bits in shots)
+
+# bits[q] is qubit q, so these are |q1 q0> = |00> and |11>.
+p00 = sim.prob_bitstring([False, False])
+p11 = sim.prob_bitstring([True, True])
+assert math.isclose(p00, 0.5, abs_tol=1e-12)
+assert math.isclose(p11, 0.5, abs_tol=1e-12)
+
+# Python reads auto-flush lazy operations and merged RZ rotations.
+accuracy = {
+    "state_exact": sim.is_state_exact(),
+    "pragmatic_drift_count": sim.pragmatic_drift_count,
+    "truncation_error": sim.truncation_error,
+    "bond_cap_hits": sim.bond_cap_hits,
+}
+assert accuracy == {
+    "state_exact": True,
+    "pragmatic_drift_count": 0,
+    "truncation_error": 0.0,
+    "bond_cap_hits": 0,
+}
+```
+
+The accuracy fields answer different questions:
+
+- `is_state_exact()` detects pending work, an unmaterialized Pauli frame, and stored-state drift from eager measurement. It does not include MPS truncation in its definition.
+- `pragmatic_drift_count` should remain zero when exact amplitudes are required after random measurements. Construct with `lazy_measure=True` to avoid that eager-measurement drift.
+- `truncation_error` estimates accumulated discarded singular-value weight.
+- `bond_cap_hits` counts SVDs at which `max_bond_dim` was binding.
+
+Python state reads automatically flush lazy operations and merged rotations. When Pauli-frame tracking is enabled, call `flush_pauli_frame_to_state()` before a read that must include the physical frame.
+
+## `Mast` quickstart
+
+`max_non_clifford` reserves one fresh ancilla for each deferred non-Clifford RZ. Exceeding it raises `PanicException`, so use compile advice to size the simulator and inspect `remaining_injections` while building a circuit. Prefer MAST for T-like gates whose injection corrections are Clifford and when the extra ancillas fit; prefer `StabMps` for direct arbitrary rotations, limited ancillary memory, or amplitude, probability, and bulk-sampling reads.
+
+```python
+import pecos_rslib_exp as exp
+
+mast = exp.Mast(2, max_non_clifford=2, seed=11)
+mast.run_gate("H", {0})
+mast.run_gate("CX", {(0, 1)})
+mast.run_gate("T", {1})
+
+assert mast.num_ancillas_used == 1
+assert mast.remaining_injections == 1
+
+# Complete all deferred injections and apply their corrections.
+mast.project_all()
+assert len(mast.projection_records()) == 1
+
+# MZ would call project_all() automatically if work were still deferred.
+outcome = mast.run_1q_gate("MZ", 0)
+assert outcome in (0, 1)
+```
+
+`flush()` is not the MAST completion operation: it materializes lazy-measurement operations and pending merged rotations, but leaves already deferred injections alone. Finish with `project_all()` or measure a data qubit with MZ.
+
+## Analyze first with `StabMpsCompile`
+
+Replay the same gates through `StabMpsCompile`, then call `recommend()` or `advise()`. Advice is heuristic. Deferred capacity counts every non-Clifford RZ, even an arbitrary-angle rotation whose eventual correction is also non-Clifford.
+
+```python
+import pecos_rslib_exp as exp
+
+analysis = exp.StabMpsCompile(20)
+analysis.run_gate("H", {0, 1})
+analysis.run_gate("CX", {(0, 2), (1, 3)})
+analysis.run_gate("T", {0, 1})
+
+recommendation = analysis.recommend()
+assert recommendation["simulator"] == "stab_mps"
+
+required = analysis.nonclifford_rz_total
+assert required == 2
+
+sufficient = analysis.advise(ancilla_budget=required)
+assert sufficient["injection"] == "deferred"
+assert sufficient["deferred_feasible"] is True
+assert sufficient["simulator"] == "mast"
+
+insufficient = analysis.advise(ancilla_budget=required - 1)
+assert insufficient["injection"] == "immediate"
+assert insufficient["deferred_feasible"] is False
+assert insufficient["warnings"]
+
+unspecified = analysis.advise()
+assert unspecified["injection"] == "deferred"
+assert unspecified["deferred_feasible"] is None
+assert unspecified["warnings"]
+```
+
+`recommend()` uses ordered thresholds: pure Clifford circuits select CH form; otherwise `n <= 14` selects a dense state vector; otherwise nullity `<= 6` selects `StabMps`; otherwise non-Clifford count `<= 40` selects `StabVec`; remaining circuits select `StabMps` with adaptive bond growth suggested. `bond_dim_bound` returns `2**nullity` and saturates at the platform maximum integer if that power overflows.
+
+## Rust quickstart
+
+Rust reads do not auto-flush. Call `flush()` before state reads when lazy measurement or merged RZ is enabled, and materialize a tracked Pauli frame separately when required.
+
+```rust
+use pecos_core::{Angle64, QubitId};
+use pecos_simulators::{ArbitraryRotationGateable, CliffordGateable};
+use pecos_stab_tn::stab_mps::StabMps;
+
+fn main() {
+    let mut sim = StabMps::builder(2)
+        .seed(7)
+        .lazy_measure(true)
+        .build();
+    sim.h(&[QubitId(0)]);
+    sim.cx(&[(QubitId(0), QubitId(1))]);
+    sim.rz(Angle64::QUARTER_TURN / 2_u64, &[QubitId(1)]);
+
+    let outcome = sim.mz(&[QubitId(0)])[0].outcome;
+    sim.flush();
+    let bits = [outcome, outcome];
+    assert!((sim.prob_bitstring(&bits) - 1.0).abs() < 1e-12);
+}
+```

exp/pecos-stab-tn/examples/advice.rs

@@ -0,0 +1,57 @@
+// Copyright 2026 The PECOS Developers
+//
+// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file
+// except in compliance with the License. You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software distributed under the
+// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
+// express or implied. See the License for the specific language governing permissions and
+// limitations under the License.
+
+//! Replay a circuit, recommend a simulator, and compare three ancilla budgets.
+
+use pecos_core::{Angle64, QubitId};
+use pecos_simulators::{ArbitraryRotationGateable, CliffordGateable};
+use pecos_stab_tn::stab_mps::compile::{ExecutionAdvice, StabMpsCompile};
+
+fn analyzed_circuit() -> StabMpsCompile {
+    let mut analysis = StabMpsCompile::new(20);
+    analysis.h(&[QubitId(0), QubitId(1)]);
+    analysis.cx(&[(QubitId(0), QubitId(2)), (QubitId(1), QubitId(3))]);
+    analysis.rz(Angle64::QUARTER_TURN / 2_u64, &[QubitId(0), QubitId(1)]);
+    analysis
+}
+
+fn print_advice(label: &str, advice: &ExecutionAdvice) {
+    println!("{label}:");
+    println!("  simulator: {:?}", advice.simulator);
+    println!("  injection: {:?}", advice.injection);
+    println!("  injectable_count: {}", advice.injectable_count);
+    println!(
+        "  deferred_ancillas_required: {}",
+        advice.deferred_ancillas_required
+    );
+    println!("  deferred_feasible: {:?}", advice.deferred_feasible);
+    println!("  warnings: {:?}", advice.warnings);
+    println!("  reason: {}", advice.reason);
+}
+
+fn main() {
+    let analysis = analyzed_circuit();
+    let recommendation = analysis.recommend();
+    println!(
+        "recommendation: {:?}: {}",
+        recommendation.kind, recommendation.reason
+    );
+
+    let required = usize::try_from(analysis.nonclifford_rz_total())
+        .expect("the example's gate count fits usize");
+    print_advice("sufficient budget", &analysis.advise(Some(required)));
+    print_advice(
+        "insufficient budget",
+        &analysis.advise(Some(required.saturating_sub(1))),
+    );
+    print_advice("unspecified budget", &analysis.advise(None));
+}

exp/pecos-stab-tn/examples/mast_projection_order.rs

@@ -0,0 +1,115 @@
+// Copyright 2026 The PECOS Developers
+//
+// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file
+// except in compliance with the License. You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software distributed under the
+// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
+// express or implied. See the License for the specific language governing permissions and
+// limitations under the License.
+
+//! Compare deferred-projection ordering for seeded random Clifford+T circuits.
+//!
+//! Circuit generation follows the `next_rng` xorshift helper and random
+//! Clifford-layer pattern in `disent_firing_rate.rs`.
+
+use pecos_core::{Angle64, QubitId};
+use pecos_simulators::{ArbitraryRotationGateable, CliffordGateable};
+use pecos_stab_tn::stab_mps::mast::{Mast, ProjectionOrder};
+use std::time::{Duration, Instant};
+
+#[derive(Clone, Copy)]
+enum CircuitGate {
+    H(usize),
+    Sz(usize),
+    Cx(usize, usize),
+    T(usize),
+}
+
+/// Same xorshift generator as `examples/disent_firing_rate.rs`.
+fn next_rng(state: &mut u64) -> u64 {
+    *state ^= *state << 13;
+    *state ^= *state >> 7;
+    *state ^= *state << 17;
+    *state
+}
+
+fn random_clifford_t_circuit(num_qubits: usize, t_count: usize, seed: u64) -> Vec<CircuitGate> {
+    let mut gates = (0..num_qubits).map(CircuitGate::H).collect::<Vec<_>>();
+    let mut rng_state = seed.wrapping_add(1);
+
+    for _ in 0..t_count {
+        for _ in 0..3 {
+            let gate_type = next_rng(&mut rng_state) % 3;
+            let q0 = (next_rng(&mut rng_state) % num_qubits as u64) as usize;
+            match gate_type {
+                0 => gates.push(CircuitGate::H(q0)),
+                1 => gates.push(CircuitGate::Sz(q0)),
+                _ => {
+                    let q1 = loop {
+                        let candidate = (next_rng(&mut rng_state) % num_qubits as u64) as usize;
+                        if candidate != q0 {
+                            break candidate;
+                        }
+                    };
+                    gates.push(CircuitGate::Cx(q0, q1));
+                }
+            }
+        }
+        let target = (next_rng(&mut rng_state) % num_qubits as u64) as usize;
+        gates.push(CircuitGate::T(target));
+    }
+    gates
+}
+
+fn run_circuit(
+    num_qubits: usize,
+    t_count: usize,
+    simulator_seed: u64,
+    gates: &[CircuitGate],
+    order: ProjectionOrder,
+) -> (usize, usize, Duration) {
+    let start = Instant::now();
+    let mut mast = Mast::with_seed(num_qubits, t_count, simulator_seed).projection_order(order);
+    let t = Angle64::QUARTER_TURN / 2u64;
+    for &gate in gates {
+        match gate {
+            CircuitGate::H(q) => mast.h(&[QubitId(q)]),
+            CircuitGate::Sz(q) => mast.sz(&[QubitId(q)]),
+            CircuitGate::Cx(control, target) => mast.cx(&[(QubitId(control), QubitId(target))]),
+            CircuitGate::T(q) => mast.rz(t, &[QubitId(q)]),
+        };
+    }
+    mast.project_all();
+    let elapsed = start.elapsed();
+    (mast.projection_peak_bond(), mast.max_bond_dim(), elapsed)
+}
+
+fn main() {
+    let circuit_seed = 0x5eed_c1ff_07d0_2026;
+    let simulator_seed = 0x5eed_5a1f_2026_0814;
+
+    println!("MAST deferred-projection order comparison");
+    println!(
+        "{:<6} {:<8} {:<10} {:>18} {:>14} {:>14}",
+        "data", "T-count", "order", "projection peak", "final bond", "wall time (s)"
+    );
+    println!("{:-<76}", "");
+
+    for num_qubits in [8usize, 16, 32] {
+        for t_count in [num_qubits, 2 * num_qubits] {
+            let gates = random_clifford_t_circuit(num_qubits, t_count, circuit_seed);
+            for order in [ProjectionOrder::Input, ProjectionOrder::MinSpan] {
+                let (projection_peak, final_bond, elapsed) =
+                    run_circuit(num_qubits, t_count, simulator_seed, &gates, order);
+                println!(
+                    "{num_qubits:<6} {t_count:<8} {order:<10?} {projection_peak:>18} \
+                     {final_bond:>14} {:>14.6}",
+                    elapsed.as_secs_f64()
+                );
+            }
+        }
+    }
+}