Symmetry-merged evolution: Trotter merging + momentum-sector CTPP

crates/ppvm-lindblad/Cargo.toml

@@ -10,6 +10,7 @@ ndarray = "0.17"
 num = "0.4.3"
 ppvm-traits     = { version = "0.1.0", path = "../ppvm-traits" }
 ppvm-pauli-word = { version = "0.1.0", path = "../ppvm-pauli-word" }
+ppvm-pauli-sum  = { version = "0.1.0", path = "../ppvm-pauli-sum" }
 rayon = "1.11"
 # Matrix-exponential action (Al-Mohy & Higham). QuSpin-rust is MIT-licensed;
 # the pinned rev is the commit that added the LICENSE file.
@@ -18,6 +19,3 @@ quspin-expm = { git = "https://github.com/QuSpin/QuSpin-rust", rev = "a0ad6c9fe2
 # we implement in `mf_expm.rs`) is not re-exported from `quspin-expm`'s root,
 # so we depend on `quspin-types` directly. Same git rev as `quspin-expm`.
 quspin-types = { git = "https://github.com/QuSpin/QuSpin-rust", rev = "a0ad6c9fe2e8063208f9ba1c6677150c993bb554" }
-
-[dev-dependencies]
-ppvm-pauli-sum = { version = "0.1.0", path = "../ppvm-pauli-sum" }

crates/ppvm-lindblad/src/basis.rs

@@ -4,12 +4,18 @@
 //! Basis-level `L*` operators: in-basis generator and off-basis leakage.
 
 use crate::Error;
+use crate::sector::Sector;
 use crate::spec::LindbladSpec;
+use crate::truncate::{cap_map_to_room, order_by_desc_mag};
 use crate::word::{Word, word_hash};
-use fxhash::{FxBuildHasher, FxHashMap};
+use fxhash::{FxBuildHasher, FxHashMap, FxHashSet};
 use num::Complex;
 use rayon::prelude::*;
 
+/// Chunk size for the leakage accumulation loops: candidates are folded
+/// into the live map (and the room-cap applied) once per chunk.
+const CHUNK_SIZE: usize = 4096;
+
 /// Build a `word → row` map for a basis assumed to contain unique Pauli
 /// words; debug-asserts the uniqueness invariant.
 pub fn build_basis_index(basis: &[Word]) -> FxHashMap<Word, u32> {
@@ -71,17 +77,7 @@ impl LindbladSpec {
         let protected_set: FxHashMap<u64, ()> =
             protected.iter().map(|w| (word_hash(w), ())).collect();
 
-        // Descending sort by |c|: process largest-magnitude contributors
-        // first so the running room-cap keeps the right entries.
-        let mut order: Vec<usize> = (0..basis.len()).collect();
-        order.sort_by(|&a, &b| {
-            coeffs[b]
-                .abs()
-                .partial_cmp(&coeffs[a].abs())
-                .unwrap_or(std::cmp::Ordering::Equal)
-        });
-
-        const CHUNK_SIZE: usize = 4096;
+        let order = order_by_desc_mag(coeffs);
         let room = max_basis.saturating_sub(basis.len());
         let n_qubits = self.n_qubits();
         let mut merged: FxHashMap<Word, f64> = FxHashMap::default();
@@ -119,21 +115,7 @@ impl LindbladSpec {
                     *merged.entry(k).or_insert(0.0) += val;
                 }
             }
-
-            // Room-cap: keep only the `room` largest-magnitude entries.
-            if merged.len() > room {
-                if room == 0 {
-                    merged.clear();
-                } else {
-                    let mut mags: Vec<f64> = merged.values().map(|v| v.abs()).collect();
-                    let k = room.min(mags.len() - 1);
-                    mags.select_nth_unstable_by(k, |a, b| {
-                        b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal)
-                    });
-                    let cutoff = mags[k];
-                    merged.retain(|_, &mut v| v.abs() >= cutoff);
-                }
-            }
+            cap_map_to_room(&mut merged, room);
         }
         // Rate-based admission: keep only candidates whose leakage rate
         // exceeds `tau_add`. `tau_add = 0` admits everything except exact
@@ -211,7 +193,6 @@ impl LindbladSpec {
         let protected_set: FxHashMap<u64, ()> =
             protected.iter().map(|w| (word_hash(w), ())).collect();
 
-        const CHUNK_SIZE: usize = 4096;
         let n_qubits = self.n_qubits();
         let mut merged: FxHashMap<Word, Complex<f64>> = FxHashMap::default();
         for chunk_start in (0..basis.len()).step_by(CHUNK_SIZE) {
@@ -253,4 +234,98 @@ impl LindbladSpec {
         }
         Ok(merged.into_iter().filter(|(_, c)| c.norm() > 0.0).collect())
     }
+
+    /// Phase-aware leakage: out-of-basis component of `L*(O_k)` where
+    /// `O_k` is the operator represented by `basis` (orbit reps) and
+    /// `coeffs` (complex coefficients in momentum `sector`).
+    ///
+    /// For each input rep `r` with coefficient `c_r`, and each output `q`
+    /// of `L*(r) = Σ_q v_q · q`:
+    /// 1. Canonicalize `q` → `(r_q, χ_k, |orbit_q|)` via
+    ///    [`Sector::canonicalize_phase`].
+    /// 2. If `r_q` NOT in `basis` and NOT in `protected`:
+    ///    `merged[r_q] += χ_k · v_q · c_r · |orbit_r| / |orbit_q|`.
+    ///
+    /// The `|orbit_r| / |orbit_q|` factor is the convention conversion
+    /// documented on [`crate::mf_expm`]'s `build_orbit_rep_cols`, so the
+    /// admitted rates are comparable to the averaged-convention
+    /// coefficients the caller holds.
+    ///
+    /// Returns `(r_q, sum)` pairs for all candidates with nonzero sum.
+    ///
+    /// This is the orbit-rep counterpart of [`Self::leakage_with_prune`],
+    /// and caps the live candidate map the same way: to the *available
+    /// room* `room = max_basis − basis.len()` (the reps we could actually
+    /// add), applied during accumulation. A large `max_basis`
+    /// (room ≥ all candidates) disables the cap — the near-exact case.
+    pub fn leakage_orbit_rep(
+        &self,
+        basis: &[Word],
+        coeffs: &[Complex<f64>],
+        protected: &[Word],
+        sector: Sector<'_>,
+        max_basis: usize,
+    ) -> Result<Vec<(Word, Complex<f64>)>, Error> {
+        if basis.len() != coeffs.len() {
+            return Err(Error::LengthMismatch {
+                what: "basis and coeffs",
+                a: basis.len(),
+                b: coeffs.len(),
+            });
+        }
+        // Membership is tested on the canonical rep `r_q`, so unlike the
+        // real path these are full-Word sets, not `word_hash` tables.
+        let in_basis: FxHashSet<&Word> = basis.iter().collect();
+        let protected_set: FxHashSet<&Word> = protected.iter().collect();
+
+        let order = order_by_desc_mag(coeffs);
+        let room = max_basis.saturating_sub(basis.len());
+        let n_qubits = self.n_qubits();
+        let mut merged: FxHashMap<Word, Complex<f64>> = FxHashMap::default();
+        for chunk_indices in order.chunks(CHUNK_SIZE) {
+            let local: Vec<Vec<(Word, Complex<f64>)>> = chunk_indices
+                .par_iter()
+                .map_init(
+                    || {
+                        (
+                            Vec::<u32>::with_capacity(n_qubits),
+                            Vec::<u32>::with_capacity(128),
+                            FxHashMap::<Word, Complex<f64>>::with_capacity_and_hasher(
+                                128,
+                                FxBuildHasher::default(),
+                            ),
+                        )
+                    },
+                    |(s1, s2, lm), &i| {
+                        let r = &basis[i];
+                        let c_r = coeffs[i];
+                        // A rep that cannot carry the sector contributes
+                        // nothing (its coefficient is identically zero).
+                        let Some(orbit_in) = sector.orbit_size(r) else {
+                            return Vec::new();
+                        };
+                        let terms = self.compute_action_terms(r, s1, s2, lm);
+                        let mut out = Vec::with_capacity(terms.len());
+                        for (q, v) in terms.iter() {
+                            let Some((r_q, phase, orbit_out)) = sector.canonicalize_phase(q) else {
+                                continue;
+                            };
+                            if !in_basis.contains(&r_q) && !protected_set.contains(&r_q) {
+                                let rate = phase * *v * c_r * (orbit_in as f64 / orbit_out as f64);
+                                out.push((r_q, rate));
+                            }
+                        }
+                        out
+                    },
+                )
+                .collect();
+            for v in local {
+                for (k, val) in v {
+                    *merged.entry(k).or_insert(Complex::new(0.0, 0.0)) += val;
+                }
+            }
+            cap_map_to_room(&mut merged, room);
+        }
+        Ok(merged.into_iter().filter(|(_, c)| c.norm() > 0.0).collect())
+    }
 }

crates/ppvm-lindblad/src/lib.rs

@@ -38,8 +38,11 @@ mod basis;
 pub mod config;
 pub mod error;
 pub(crate) mod expm;
+mod scalar;
+pub mod sector;
 mod spec;
 mod step;
+mod truncate;
 mod word;
 
 /// Matrix-free / quspin-expm-backed `exp(dt·L*)·b` engine. See module docs.
@@ -48,6 +51,7 @@ pub(crate) mod mf_expm;
 pub use basis::build_basis_index;
 pub use config::PcStepConfig;
 pub use error::Error;
+pub use sector::{Sector, canonicalize_basis_to_rep};
 pub use spec::{JumpInput, LindbladSpec};
 pub use step::PcStepTimings;
 pub use word::{MAX_QUBITS, Word, codes_from_word, parse_pauli_string, word_from_codes};