From e82570fef10b15487bb36a72e2b81a6198998505 Mon Sep 17 00:00:00 2001 From: alexschuckert Date: Thu, 16 Jul 2026 13:18:25 +0100 Subject: [PATCH 01/13] =?UTF-8?q?feat:=20symmetry-merged=20evolution=20?= =?UTF-8?q?=E2=80=94=20Trotter=20merging=20+=20momentum-sector=20CTPP?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two consumers of the translation-symmetry primitive: Trotter path: PauliSum.symmetry_merge (k=0) and PauliSum.momentum_merge (k≠0 carried as a real pair, character-weighted fold, |G|-rescaled to the summing projector so merging after every step is idempotent), with TranslationGroup and canonicalize_basis_arr{,_complex} / check_momentum_sector_arr exposed to Python. Inputs validated at the boundary (ValueError, not panics). CTPP path: pc_step_orbit_rep — per-step evolution entirely in orbit-representative form (complex coefficients, phase-aware action, cached-CSC expm; basis ~|G|× smaller than full-basis evolution, persisting through every step), with the same PcStepConfig truncation policy incl. displacement admission. Split 3/4 of the CTPP work; full history on branch continuous-time-pauli-propagation. Co-Authored-By: Claude Fable 5 --- crates/ppvm-lindblad/Cargo.toml | 4 +- crates/ppvm-lindblad/src/lib.rs | 5 + crates/ppvm-lindblad/src/orbit_rep.rs | 445 +++++++++++++++++++++ crates/ppvm-lindblad/src/tests.rs | 102 +++++ crates/ppvm-python-native/src/interface.rs | 139 +++++++ crates/ppvm-python-native/src/lib.rs | 11 + crates/ppvm-python-native/src/lindblad.rs | 118 +++++- crates/ppvm-python-native/src/symmetry.rs | 328 +++++++++++++++ ppvm-python/src/ppvm/lindblad.py | 60 ++- ppvm-python/src/ppvm/paulisum.py | 44 ++ ppvm-python/test/test_momentum_merge.py | 199 +++++++++ 11 files changed, 1450 insertions(+), 5 deletions(-) create mode 100644 crates/ppvm-lindblad/src/orbit_rep.rs create mode 100644 crates/ppvm-python-native/src/symmetry.rs create mode 100644 ppvm-python/test/test_momentum_merge.py diff --git a/crates/ppvm-lindblad/Cargo.toml b/crates/ppvm-lindblad/Cargo.toml index 69f425341..900388d36 100644 --- a/crates/ppvm-lindblad/Cargo.toml +++ b/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" } diff --git a/crates/ppvm-lindblad/src/lib.rs b/crates/ppvm-lindblad/src/lib.rs index da5414e66..e99bbb5d0 100644 --- a/crates/ppvm-lindblad/src/lib.rs +++ b/crates/ppvm-lindblad/src/lib.rs @@ -45,6 +45,10 @@ mod word; /// Matrix-free / quspin-expm-backed `exp(dt·L*)·b` engine. See module docs. pub(crate) mod mf_expm; +/// Per-step orbit-rep evolution under translation symmetry, with a +/// phase-aware complex action. See module docs. +pub mod orbit_rep; + pub use basis::build_basis_index; pub use config::PcStepConfig; pub use error::Error; @@ -54,3 +58,4 @@ pub use word::{MAX_QUBITS, Word, codes_from_word, parse_pauli_string, word_from_ #[cfg(test)] mod tests; + diff --git a/crates/ppvm-lindblad/src/orbit_rep.rs b/crates/ppvm-lindblad/src/orbit_rep.rs new file mode 100644 index 000000000..7bf7a27c6 --- /dev/null +++ b/crates/ppvm-lindblad/src/orbit_rep.rs @@ -0,0 +1,445 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Per-step orbit-representative evolution under translation symmetry. +//! +//! The state lives entirely in **orbit-rep form** throughout: `basis` +//! contains only canonical translation-orbit representatives, and +//! `coeffs` are complex (one per rep). The dynamics `L*` is computed +//! with **phase-aware action** — for each output Pauli `q`, we +//! canonicalize `q` to its orbit rep `r_q` with shift counter `cnt_q`, +//! and accumulate `χ_k(g_{cnt_q}) · v · c_r` (where `v` is the matrix +//! element of `L*` between input rep `r` and output `q`). +//! +//! The orbit-rep basis is ~|G|× smaller than the full-basis +//! representation, throughout the entire evolution. +//! +//! The phase-aware action is genuinely **complex** (because of the +//! `χ_k(g)` phase factors). Rather than materialise a sparse matrix, the +//! per-column action — for each input rep, the list of `(row, χ_k·v)` +//! pairs for the in-basis outputs — is computed **once per expm call** +//! (via [`build_orbit_rep_cols`]) and then reused, CSC-style, across +//! every Krylov–Taylor matvec driving the external `quspin-expm` engine. +//! +//! ## Limitations +//! +//! - Caller is responsible for ensuring the input basis is in orbit-rep +//! form (i.e. each entry is the canonical representative of its +//! translation orbit). Use [`canonicalize_basis_to_rep`] if needed. +//! - The momentum sector `k_modes` is fixed for the duration of one +//! pc_step call. To compute a full site-resolved profile, call +//! `pc_step_orbit_rep` once per momentum mode and inverse-Fourier +//! the results. + +use crate::mf_expm::CscOp; +use crate::{Error, LindbladSpec}; +use fxhash::{FxBuildHasher, FxHashMap}; +use num::Complex; +use ppvm_pauli_sum::symmetry::TranslationGroup; +use quspin_expm::ExpmOp; +use rayon::prelude::*; + +// Word type re-exported from lib.rs. +use crate::Word; + +/// Replace each entry of `basis` with its canonical orbit +/// representative under `group`. Pure rewrite; coefficients are +/// untouched. Useful to enforce the orbit-rep invariant before calling +/// [`pc_step_orbit_rep`]. +/// +/// Does NOT deduplicate — if multiple input entries collapse to the +/// same rep, both are kept (caller should run a merge afterwards). +pub fn canonicalize_basis_to_rep(basis: &mut [Word], group: &TranslationGroup) { + for w in basis.iter_mut() { + *w = group.canonicalize(w); + } +} + +/// Build the per-column phase-aware action of the in-basis-restricted +/// orbit-rep generator `M` at momentum sector `k_modes`. +/// +/// Returns, for each input rep `basis[c]` (column `c`), the list of +/// `(row, χ_k(g_{cnt_q}) · v_q)` pairs for every action output Pauli `q` +/// of `L*(basis[c])` whose orbit rep `r_q` is in `basis` at index `row`. +/// Outputs not in `basis` are dropped. This is the expensive part of the +/// orbit-rep dynamics (`compute_action_terms`, `canonicalize_with_shift`, +/// `character`); it is computed once and reused by the CSC-style matvec +/// in [`CscOp`]. +pub(crate) fn build_orbit_rep_cols( + spec: &LindbladSpec, + basis: &[Word], + index: &FxHashMap, + group: &TranslationGroup, + k_modes: &[i32], +) -> Vec)>> { + basis + .par_iter() + .map_init( + || { + ( + Vec::::with_capacity(spec.n_qubits()), + Vec::::with_capacity(128), + FxHashMap::>::with_capacity_and_hasher( + 128, + FxBuildHasher::default(), + ), + ) + }, + |(s1, s2, lm), r| { + let terms = spec.compute_action_terms(r, s1, s2, lm); + let mut out = Vec::with_capacity(terms.len()); + for (q, v) in terms.iter() { + let (r_q, cnt_q) = group.canonicalize_with_shift(q); + if let Some(&row) = index.get(&r_q) { + let phase = group.character(k_modes, &cnt_q); + out.push((row, phase * *v)); + } + } + out + }, + ) + .collect() +} + +/// Compute `exp(dt · M) · coeffs` for the in-basis-restricted orbit-rep +/// generator `M` at momentum sector `k_modes`, via `quspin-expm`. Returns +/// a fresh `Vec>` of length `basis.len()`. +/// +/// The expensive phase-aware action is computed ONCE here (via +/// [`build_orbit_rep_cols`]) and reused, CSC-style, across every Krylov– +/// Taylor matvec (see [`CscOp`]). One pass over the cached columns +/// extracts the diagonal shift `μ = tr(M)/n` and a valid upper bound on +/// the column 1-norm of `M − μ·I`; from `‖dt·(M−μI)‖₁` we pick the Taylor +/// partition `(m*, s)` and hand everything to +/// [`quspin_expm::ExpmOp::from_parts`] (mirroring +/// [`crate::mf_expm::expm_apply_mf`]). +pub(crate) fn expm_apply_orbit_rep_cached( + spec: &LindbladSpec, + basis: &[Word], + group: &TranslationGroup, + k_modes: &[i32], + dt: f64, + coeffs: &[Complex], +) -> Vec> { + let n = basis.len(); + if n == 0 { + return Vec::new(); + } + + let index = crate::build_basis_index(basis); + let cols = build_orbit_rep_cols(spec, basis, &index, group, k_modes); + + // One pass for the per-column `(raw, diag)` used by the `μ`/1-norm + // selection: `raw = Σ|val|` (upper bound on the absolute column sum), + // `diag = M[c,c]`. From these: `trace = Σ diag`, `μ = trace/n`, and an + // upper bound on the column 1-norm of `M − μ·I`: `raw − |diag| + |diag − μ|`. + let per_col: Vec<(f64, Complex)> = cols + .par_iter() + .enumerate() + .map(|(c, col)| { + let mut raw = 0.0_f64; + let mut diag = Complex::new(0.0, 0.0); + for &(row, val) in col.iter() { + raw += val.norm(); + if row as usize == c { + diag += val; + } + } + (raw, diag) + }) + .collect(); + + let trace: Complex = per_col.iter().map(|(_, d)| *d).sum(); + let mu = trace / n as f64; + let onenorm = per_col + .iter() + .map(|(raw, diag)| raw - diag.norm() + (diag - mu).norm()) + .fold(0.0_f64, f64::max); + + let (m_star, s) = crate::expm::select_ms(dt.abs() * onenorm); + + let mut v = coeffs.to_vec(); + let op = CscOp { cols: &cols, dim: n }; + let e = ExpmOp::from_parts( + op, + Complex::new(dt, 0.0), + mu, + s as usize, + m_star as usize, + 1e-12_f64, + ); + e.apply(ndarray::ArrayViewMut1::from(v.as_mut_slice())) + .expect("expm apply (orbit-rep cached)"); + v +} + +/// 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 `k_modes`). +/// +/// 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, cnt_q)`. +/// 2. If `r_q` NOT in `basis` and NOT in `protected`: +/// `merged[r_q] += χ_k(g_{cnt_q}) · v_q · c_r`. +/// +/// Returns `(r_q, sum)` pairs for all candidates with nonzero sum. +/// +/// The live candidate map is capped to the *available room* +/// `room = max_basis − basis.len()` (the reps we could actually add), +/// applied during accumulation: input reps are processed in descending +/// `|c|` order and after each chunk only the `room` largest-`|sum|` +/// candidates are kept. A large `max_basis` (room ≥ all candidates) +/// disables the cap — the near-exact case. +#[allow(clippy::too_many_arguments)] +pub fn leakage_orbit_rep( + spec: &LindbladSpec, + basis: &[Word], + coeffs: &[Complex], + protected: &[Word], + group: &TranslationGroup, + k_modes: &[i32], + max_basis: usize, +) -> Result)>, Error> { + if basis.len() != coeffs.len() { + return Err(Error::LengthMismatch { + what: "basis and coeffs", + a: basis.len(), + b: coeffs.len(), + }); + } + let in_basis: FxHashMap<&Word, ()> = basis.iter().map(|w| (w, ())).collect(); + let protected_set: FxHashMap<&Word, ()> = protected.iter().map(|w| (w, ())).collect(); + + // Descending sort by |c|: process largest-magnitude contributors first + // so the running room-cap keeps the right entries. + let mut order: Vec = (0..basis.len()).collect(); + order.sort_by(|&a, &b| { + coeffs[b] + .norm() + .partial_cmp(&coeffs[a].norm()) + .unwrap_or(std::cmp::Ordering::Equal) + }); + + const CHUNK_SIZE: usize = 4096; + let room = max_basis.saturating_sub(basis.len()); + let mut merged: FxHashMap> = FxHashMap::default(); + for chunk_indices in order.chunks(CHUNK_SIZE) { + let local: Vec)>> = chunk_indices + .par_iter() + .map_init( + || { + ( + Vec::::with_capacity(spec.n_qubits()), + Vec::::with_capacity(128), + FxHashMap::>::with_capacity_and_hasher( + 128, + FxBuildHasher::default(), + ), + ) + }, + |(s1, s2, lm), &i| { + let r = &basis[i]; + let c_r = coeffs[i]; + let terms = spec.compute_action_terms(r, s1, s2, lm); + let mut out = Vec::with_capacity(terms.len()); + for (q, v) in terms.iter() { + let (r_q, cnt_q) = group.canonicalize_with_shift(q); + if !in_basis.contains_key(&r_q) && !protected_set.contains_key(&r_q) { + let phase = group.character(k_modes, &cnt_q); + out.push((r_q, phase * *v * c_r)); + } + } + out + }, + ) + .collect(); + for v in local { + for (k, val) in v { + *merged.entry(k).or_insert(Complex::new(0.0, 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 = merged.values().map(|v| v.norm()).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.norm() >= cutoff); + } + } + } + Ok(merged.into_iter().filter(|(_, c)| c.norm() > 0.0).collect()) +} + +/// Per-step orbit-rep predictor-corrector evolution. +/// +/// All state lives in orbit-rep form throughout. Each pc step does: +/// 1. Phase-aware leakage from `(basis, coeffs)`; append the largest +/// leakage reps, up to the admission room. +/// 2. Predictor: cached-action expm ([`expm_apply_orbit_rep_cached`]). +/// 3. Phase-aware leakage from the predicted state; append further reps. +/// 4. Corrector: cached-action expm from the pre-step coefficients on +/// the doubly-enlarged basis. +/// 5. Prune `|c| < drop_tol`, then trim to the top-`max_basis` reps by +/// `|c|`; protected reps never dropped. +/// +/// `max_basis` is a hard rank cap on the live orbit-rep basis: enrichment +/// adds at most `max_basis − basis.len()` of the largest leakage reps, the +/// leakage map is capped to the same room, and the post-step basis is +/// trimmed to the top-`max_basis` by `|c|`. Pass a large value (e.g. +/// `usize::MAX`) for the near-exact, uncapped case. `drop_tol` additionally +/// prunes by magnitude. +/// +/// `basis` is assumed to contain only canonical orbit representatives. +/// If not, [`canonicalize_basis_to_rep`] should be called first. +#[allow(clippy::too_many_arguments)] +pub fn pc_step_orbit_rep( + spec: &LindbladSpec, + basis: &mut Vec, + coeffs: &mut Vec>, + dt: f64, + protected: &[Word], + group: &TranslationGroup, + k_modes: &[i32], + cfg: &crate::PcStepConfig, +) -> Result<(), Error> { + let crate::PcStepConfig { max_basis, admit_basis, drop_tol, tau_add, .. } = *cfg; + // Admission bound, mirroring the real-space `pc_step`: enrichment may + // grow the live basis to `admit` >= `max_basis`; the final + // `cap_basis_complex` keeps the top-`max_basis` reps by evolved |coeff| + // over the whole union (rank displacement). With `admit_basis = None` + // admission is bounded by `max_basis` itself and membership turnover + // requires `drop_tol > 0`. + let admit = admit_basis.unwrap_or(max_basis).max(max_basis); + let tau_add = tau_add.unwrap_or(0.0); + // 1. First-hop phase-aware leakage. + let mut leak = leakage_orbit_rep(spec, basis, coeffs, protected, group, k_modes, admit)?; + if tau_add > 0.0 { + leak.retain(|(_, c)| c.norm() > tau_add); + } + add_leakage_capped_complex(basis, coeffs, leak, admit); + // 2. Predictor: cached-action expm (the phase-aware action is built + // once via `build_orbit_rep_cols`). + let coeffs_predict = expm_apply_orbit_rep_cached(spec, basis, group, k_modes, dt, coeffs); + // 3. Second-hop leakage from predicted state. + let mut leak2 = + leakage_orbit_rep(spec, basis, &coeffs_predict, protected, group, k_modes, admit)?; + drop(coeffs_predict); + if tau_add > 0.0 { + leak2.retain(|(_, c)| c.norm() > tau_add); + } + add_leakage_capped_complex(basis, coeffs, leak2, admit); + // 4. Corrector: cache-the-action expm from pre-step state (basis grew). + *coeffs = expm_apply_orbit_rep_cached(spec, basis, group, k_modes, dt, coeffs); + // 5. Prune by magnitude, then rank-cap to max_basis. + if drop_tol > 0.0 { + prune_basis_complex_local(basis, coeffs, drop_tol, protected); + } + cap_basis_complex(basis, coeffs, max_basis, protected); + Ok(()) +} + +/// Complex analogue of `crate::add_leakage_capped`: add the largest leakage +/// reps to the basis, up to the available room `room = max_basis − +/// basis.len()`, so the in-step orbit-rep basis never exceeds `max_basis`. +/// New reps get coefficient 0; the surrounding expm fills them. No +/// magnitude filter — the top-`room` by `|leakage|` are added. +fn add_leakage_capped_complex( + basis: &mut Vec, + coeffs: &mut Vec>, + mut leak: Vec<(Word, Complex)>, + max_basis: usize, +) { + let room = max_basis.saturating_sub(basis.len()); + if leak.len() > room { + if room > 0 { + leak.select_nth_unstable_by(room - 1, |a, b| { + b.1.norm() + .partial_cmp(&a.1.norm()) + .unwrap_or(std::cmp::Ordering::Equal) + }); + } + leak.truncate(room); + } + for (w, _) in leak { + basis.push(w); + coeffs.push(Complex::new(0.0, 0.0)); + } +} + +/// Complex analogue of `crate::cap_basis`: keep only the `max_basis` +/// largest-`|c|` reps (protected reps always kept), dropping the rest. +/// A `max_basis` large enough to cover the whole basis is a no-op. +fn cap_basis_complex( + basis: &mut Vec, + coeffs: &mut Vec>, + max_basis: usize, + protected: &[Word], +) { + if basis.len() <= max_basis { + return; + } + let protected_set: fxhash::FxHashSet<&Word> = protected.iter().collect(); + let n_prot = basis.iter().filter(|w| protected_set.contains(w)).count(); + let slots = max_basis.saturating_sub(n_prot); + let mut mags: Vec = basis + .iter() + .zip(coeffs.iter()) + .filter(|(w, _)| !protected_set.contains(w)) + .map(|(_, c)| c.norm()) + .collect(); + let cutoff = if slots == 0 { + f64::INFINITY + } else if slots >= mags.len() { + return; + } else { + let k = slots - 1; + mags.select_nth_unstable_by(k, |a, b| b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal)); + mags[k] + }; + let mut write = 0; + for read in 0..basis.len() { + if protected_set.contains(&basis[read]) || coeffs[read].norm() >= cutoff { + if write != read { + basis.swap(write, read); + coeffs.swap(write, read); + } + write += 1; + } + } + basis.truncate(write); + coeffs.truncate(write); +} + +/// Complex analogue of `crate::prune_basis`: drop reps with `|c| < +/// drop_tol`, never dropping `protected` reps. No-op when `drop_tol <= 0`. +fn prune_basis_complex_local( + basis: &mut Vec, + coeffs: &mut Vec>, + drop_tol: f64, + protected: &[Word], +) { + if drop_tol <= 0.0 { + return; + } + let protected_set: fxhash::FxHashSet<&Word> = protected.iter().collect(); + let mut write = 0; + for read in 0..basis.len() { + if coeffs[read].norm() >= drop_tol || protected_set.contains(&basis[read]) { + if write != read { + basis.swap(write, read); + coeffs.swap(write, read); + } + write += 1; + } + } + basis.truncate(write); + coeffs.truncate(write); +} diff --git a/crates/ppvm-lindblad/src/tests.rs b/crates/ppvm-lindblad/src/tests.rs index 4217b3e13..3ad4899d5 100644 --- a/crates/ppvm-lindblad/src/tests.rs +++ b/crates/ppvm-lindblad/src/tests.rs @@ -94,6 +94,108 @@ fn word_codec_roundtrip() { assert_eq!(out.as_slice(), &codes); } +/// Per-step orbit-rep evolution gives the SAME final orbit-rep +/// state as full-basis complex evolution followed by a single +/// projection at the end. Validates that the phase-aware complex +/// action machinery is consistent with the full-basis reference. +#[test] +fn pc_step_orbit_rep_matches_full_basis_projection() { + use std::f64::consts::PI; + + use ppvm_pauli_sum::symmetry::canonicalize_pauli_sum_complex; + let n = 4usize; + let dt = 0.01f64; + let n_steps = 3usize; + let mut h_terms: Vec<(String, f64)> = Vec::new(); + for j in 0..n { + let nxt = (j + 1) % n; + for op in ["X", "Y"] { + let mut s = vec!['I'; n]; + s[j] = op.chars().next().unwrap(); + s[nxt] = op.chars().next().unwrap(); + h_terms.push((s.into_iter().collect(), 1.0)); + } + } + let spec = LindbladSpec::new(n, &h_terms, &[]).unwrap(); + let group = ppvm_pauli_sum::symmetry::TranslationGroup::chain_1d(n); + let k_mode: i32 = 1; + let k = vec![k_mode]; + + // Build the k=1 eigenstate in FULL basis form. + let basis_full: Vec = (0..n) + .map(|j| { + let mut s = vec!['I'; n]; + s[j] = 'Z'; + let (w, _) = parse_pauli_string(&s.into_iter().collect::(), n).unwrap(); + w + }) + .collect(); + let coeffs_full: Vec> = (0..n as i32) + .map(|a| Complex::from_polar(1.0, -2.0 * PI * (k_mode as f64) * (a as f64) / (n as f64))) + .collect(); + + // ----- Full-basis path ----- + let mut bf = basis_full.clone(); + let mut cf = coeffs_full.clone(); + let protected: Vec = Vec::new(); + for _ in 0..n_steps { + // Full enrichment (tau_add = 0.0 adds every leakage string): + // for a momentum eigenstate the leakage is pure-sector, so the + // full-basis and orbit-rep paths build corresponding bases and + // the projection theorem gives an exact match. The orbit-rep + // side uses a large max_basis so its rank cap never binds. + pc_step_complex_full(&spec, &mut bf, &mut cf, dt); + } + // Project at the end. + canonicalize_pauli_sum_complex(&mut bf, &mut cf, &group, &k); + + // ----- Orbit-rep path ----- + // Initial orbit-rep form: project the full-basis input. + let mut br = basis_full.clone(); + let mut cr = coeffs_full.clone(); + canonicalize_pauli_sum_complex(&mut br, &mut cr, &group, &k); + // Evolve in orbit-rep form (max_basis large ⇒ full enrichment). + for _ in 0..n_steps { + orbit_rep::pc_step_orbit_rep( + &spec, + &mut br, + &mut cr, + dt, + &protected, + &group, + &k, + &PcStepConfig { + max_basis: 10_000_000, + ..Default::default() + }, + ) + .unwrap(); + } + + // Compare. + let mf: FxHashMap> = bf.into_iter().zip(cf).collect(); + let mr: FxHashMap> = br.into_iter().zip(cr).collect(); + assert_eq!( + mf.len(), + mr.len(), + "orbit-rep ({}) and full-basis-projected ({}) basis sizes differ", + mr.len(), + mf.len() + ); + let mut max_diff = 0.0_f64; + for (w, cm) in &mr { + let cf_val = mf + .get(w) + .copied() + .unwrap_or_else(|| panic!("rep {:?} in orbit-rep but not in full-basis", w)); + max_diff = max_diff.max((cm - cf_val).norm()); + } + assert!( + max_diff < 1e-9, + "orbit-rep diverged from full-basis: max |Δc| = {max_diff:e}" + ); +} + /// The full-space complex step at momentum k=0 must reproduce the real /// pc_step on the same trajectory exactly. #[test] diff --git a/crates/ppvm-python-native/src/interface.rs b/crates/ppvm-python-native/src/interface.rs index 44a47b83d..ec690ae80 100644 --- a/crates/ppvm-python-native/src/interface.rs +++ b/crates/ppvm-python-native/src/interface.rs @@ -48,6 +48,139 @@ macro_rules! create_interface_loss_methods { }; } +macro_rules! create_interface_symmetry_methods { + // Skip loss variants: LossyPauliWord canonicalization would need + // simultaneous permutation of the loss bitmap, which we don't + // implement here. + ($name: ident, $type: ident, true) => {}; + ($name: ident, $type: ident, false) => { + #[pymethods] + impl $name { + /// Symmetry-merge this PauliSum in place: replace every + /// Pauli word by its canonical orbit representative under + /// `group`, accumulating coefficients on collision. Reduces + /// entry count by up to `|group|×` for translation-invariant + /// operators. + /// + /// See `ppvm._core.TranslationGroup` for constructors + /// (`chain_1d`, `torus_2d`, `torus_3d`, `ladder`). + /// + /// Plain real-coefficient merge (the `k=0` symmetry sector). + /// For non-trivial momentum sectors use `momentum_merge`. + pub fn symmetry_merge( + &mut self, + group: &crate::symmetry::TranslationGroup, + ) -> pyo3::PyResult<()> { + if self.inner.n_qubits() != group.core().n_qubits() { + return Err(pyo3::exceptions::PyValueError::new_err(format!( + "PauliSum has {} qubits but the TranslationGroup acts on {}", + self.inner.n_qubits(), + group.core().n_qubits(), + ))); + } + ppvm_pauli_sum::symmetry::symmetry_merge_pauli_sum( + &mut self.inner, + group.core(), + ); + Ok(()) + } + + /// Phase-aware (momentum-sector) merge for a complex operator + /// carried as a *real pair*: `self` is the real part, `other` + /// the imaginary part of `O = self + i·other`. Both are + /// overwritten in place with the orbit-representative form + /// projected onto momentum sector `momentum` (one integer mode + /// per group generator; `[0,…]` is the trivial sector and + /// reduces to `symmetry_merge`). This generalizes + /// `symmetry_merge` to k != 0 while keeping real coefficients on + /// the Python side — the only place complex arithmetic appears + /// is the internal character-weighted fold, reusing the tested + /// `canonicalize_pauli_sum_complex`. + /// + /// `self` and `other` must be distinct objects with identical + /// qubit count. After a translation-covariant gate layer this + /// is exact; under a generic Trotter step it carries the same + /// O(dt^{p+1}) equivariance error as the k=0 merge. + #[pyo3(signature = (other, group, momentum))] + pub fn momentum_merge( + &mut self, + mut other: pyo3::PyRefMut<'_, Self>, + group: &crate::symmetry::TranslationGroup, + momentum: Vec, + ) -> pyo3::PyResult<()> { + let n_g = group.core().n_qubits(); + for (label, n) in [ + ("self", self.inner.n_qubits()), + ("other", other.inner.n_qubits()), + ] { + if n != n_g { + return Err(pyo3::exceptions::PyValueError::new_err(format!( + "{label} PauliSum has {n} qubits but the \ + TranslationGroup acts on {n_g}", + ))); + } + } + if momentum.len() != group.core().n_generators() { + return Err(pyo3::exceptions::PyValueError::new_err(format!( + "momentum has {} entries but the group has {} generators", + momentum.len(), + group.core().n_generators(), + ))); + } + // Gather both real components into word -> (re + i·im). + let mut combined: std::collections::HashMap< + <$type as Config>::PauliWordType, + num::Complex, + > = std::collections::HashMap::new(); + for (w, v) in self.inner.data().iter() { + combined + .entry(w.clone()) + .or_insert(num::Complex::new(0.0, 0.0)) + .re += *v; + } + for (w, v) in other.inner.data().iter() { + combined + .entry(w.clone()) + .or_insert(num::Complex::new(0.0, 0.0)) + .im += *v; + } + let mut basis = Vec::with_capacity(combined.len()); + let mut coeffs = Vec::with_capacity(combined.len()); + for (w, c) in combined { + basis.push(w); + coeffs.push(c); + } + // Character-weighted fold onto orbit reps. + // `canonicalize_pauli_sum_complex` carries a 1/|G| prefactor; + // we rescale by |G| so the merge is the *summing* projector + // (like `symmetry_merge`): idempotent on already-merged input, + // hence stable under merging after every Trotter step. + ppvm_pauli_sum::symmetry::canonicalize_pauli_sum_complex( + &mut basis, + &mut coeffs, + group.core(), + &momentum, + ); + let scale = group.core().order() as f64; + // Write the real/imag parts back into the two sums. + self.inner.data_mut().clear(); + other.inner.data_mut().clear(); + for (w, c) in basis.into_iter().zip(coeffs.into_iter()) { + let re = c.re * scale; + let im = c.im * scale; + if re != 0.0 { + self.inner += (w.clone(), re); + } + if im != 0.0 { + other.inner += (w, im); + } + } + Ok(()) + } + } + }; +} + macro_rules! create_strategy { (false, $min_abs_coeff:ident, $max_pauli_weight:ident, $_max_loss_weight:ident) => { CombinedStrategy( @@ -403,6 +536,12 @@ macro_rules! create_interface { } } + // `symmetry_merge` only makes sense on non-loss variants — the + // canonicalization permutes qubit positions and the loss + // bitmap would need a parallel permutation that we don't + // attempt here. + create_interface_symmetry_methods!($name, $type, $loss); + create_interface_loss_methods!($name, $type, $loss); }; } diff --git a/crates/ppvm-python-native/src/lib.rs b/crates/ppvm-python-native/src/lib.rs index 13f4c22ec..efb190104 100644 --- a/crates/ppvm-python-native/src/lib.rs +++ b/crates/ppvm-python-native/src/lib.rs @@ -16,6 +16,7 @@ pub mod interface_tableau; pub mod interface_tableau_sum; pub mod lindblad; pub mod stim_program; +pub mod symmetry; pub(crate) fn flat_pairs(targets: &[usize]) -> PyResult> { if !targets.len().is_multiple_of(2) { @@ -310,4 +311,14 @@ pub mod _core { // Lindbladian time evolution #[pymodule_export] pub use crate::lindblad::LindbladSpec; + + // Symmetry merging + #[pymodule_export] + pub use crate::symmetry::TranslationGroup; + #[pymodule_export] + pub use crate::symmetry::canonicalize_basis_arr; + #[pymodule_export] + pub use crate::symmetry::canonicalize_basis_arr_complex; + #[pymodule_export] + pub use crate::symmetry::check_momentum_sector_arr; } diff --git a/crates/ppvm-python-native/src/lindblad.rs b/crates/ppvm-python-native/src/lindblad.rs index 25d4e0b0d..07d0e739f 100644 --- a/crates/ppvm-python-native/src/lindblad.rs +++ b/crates/ppvm-python-native/src/lindblad.rs @@ -13,7 +13,9 @@ use std::collections::HashMap; use num::Complex; -use numpy::{IntoPyArray, PyArray1, PyArray2, PyArrayMethods, PyReadonlyArray1, PyReadonlyArray2}; +use numpy::{ + Complex64, IntoPyArray, PyArray1, PyArray2, PyArrayMethods, PyReadonlyArray1, PyReadonlyArray2, +}; use ppvm_lindblad::{JumpInput, LindbladSpec as CoreSpec, Word, codes_from_word, word_from_codes}; use pyo3::{exceptions::PyValueError, prelude::*}; @@ -341,6 +343,120 @@ impl LindbladSpec { Ok((map, d)) } + /// Per-step orbit-rep predictor-corrector evolution under + /// translation symmetry. State lives entirely in **orbit-rep form**: + /// basis contains only canonical orbit representatives, coefficients + /// are complex. The action is phase-aware: output Paulis canonicalize + /// to their orbit rep with momentum-character weight. + /// + /// Per-step memory benefit: basis is ~|group|× smaller than the + /// full-basis representation, and the reduction persists through + /// every step. + /// + /// **Pre-condition**: every row of `basis` must be the canonical + /// orbit representative of its translation orbit under `group`. + /// Pass `canonicalize_first=True` to enforce this on entry (rewrites + /// each basis row to its canonical rep; coefficients unchanged). + /// Default `False` — the caller is trusted. + /// + /// `max_basis` is a hard rank cap on the live orbit-rep basis: + /// enrichment adds at most `max_basis − basis.len()` of the largest + /// leakage reps and the post-step basis is trimmed to the top-`max_basis` + /// by `|c|` (protected reps always kept). Pass a large value for the + /// near-exact case. `drop_tol` additionally prunes by magnitude. + #[pyo3(signature = ( + basis, coeffs, dt, max_basis, + group, momentum, + drop_tol = 0.0, + protected = None, + canonicalize_first = false, + admit_basis = None, + tau_add = None, + ))] + #[allow(clippy::too_many_arguments)] + fn pc_step_orbit_rep<'py>( + &self, + py: Python<'py>, + basis: PyReadonlyArray2<'py, u8>, + coeffs: PyReadonlyArray1<'py, Complex64>, + dt: f64, + max_basis: usize, + group: &crate::symmetry::TranslationGroup, + momentum: PyReadonlyArray1<'py, i32>, + drop_tol: f64, + protected: Option>, + canonicalize_first: bool, + admit_basis: Option, + tau_add: Option, + ) -> PyResult<(Bound<'py, PyArray2>, Bound<'py, PyArray1>)> { + use num::Complex; + use ppvm_lindblad::orbit_rep; + + let n_q = self.inner.n_qubits(); + let basis_view = basis.as_array(); + let mut basis_words = decode_basis(&basis_view, n_q)?; + let coeffs_slice = coeffs.as_slice()?; + if coeffs_slice.len() != basis_words.len() { + return Err(PyValueError::new_err(format!( + "coeffs has length {} but basis has {} rows", + coeffs_slice.len(), + basis_words.len() + ))); + } + let mut coeffs_vec: Vec> = coeffs_slice + .iter() + .map(|c| Complex::new(c.re, c.im)) + .collect(); + let protected_words: Vec = if let Some(ref p) = protected { + decode_basis(&p.as_array(), n_q)? + } else { + Vec::new() + }; + let k_slice = momentum.as_slice()?; + if k_slice.len() != group.core().n_generators() { + return Err(PyValueError::new_err(format!( + "momentum has {} entries but group has {} generators", + k_slice.len(), + group.core().n_generators() + ))); + } + if canonicalize_first { + orbit_rep::canonicalize_basis_to_rep(&mut basis_words, group.core()); + } + orbit_rep::pc_step_orbit_rep( + &self.inner, + &mut basis_words, + &mut coeffs_vec, + dt, + &protected_words, + group.core(), + k_slice, + &ppvm_lindblad::PcStepConfig { + max_basis, + admit_basis, + drop_tol, + tau_add, + num_threads: None, + }, + ) + .map_err(map_err)?; + + let m = basis_words.len(); + let mut out_basis = vec![0u8; m * n_q]; + for (i, w) in basis_words.iter().enumerate() { + codes_from_word(w, &mut out_basis[i * n_q..(i + 1) * n_q]); + } + let out_coeffs: Vec = coeffs_vec + .iter() + .map(|c| Complex64::new(c.re, c.im)) + .collect(); + let basis_arr = out_basis + .into_pyarray(py) + .reshape([m, n_q]) + .map_err(|e| PyValueError::new_err(format!("reshape failed: {e}")))?; + Ok((basis_arr, out_coeffs.into_pyarray(py))) + } + /// Sparse generator matrix in COO form: `(rows, cols, vals)`. fn generator<'py>( &self, diff --git a/crates/ppvm-python-native/src/symmetry.rs b/crates/ppvm-python-native/src/symmetry.rs new file mode 100644 index 000000000..d3c69b6e9 --- /dev/null +++ b/crates/ppvm-python-native/src/symmetry.rs @@ -0,0 +1,328 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Python bindings for the symmetry-merging primitive. +//! +//! Exposes: +//! - [`TranslationGroup`] PyO3 class with constructors for 1D, 2D, 3D +//! tori and multi-leg ladders, plus a generic generator-list path. +//! - [`canonicalize_basis_arr`] / [`canonicalize_basis_arr_complex`] free +//! functions that merge the numpy `(basis_arr, coeffs)` representation +//! used by `Lindbladian.pc_step_arr`. + +use num::Complex; +use numpy::{ + Complex64, IntoPyArray, PyArray1, PyArray2, PyArrayMethods, PyReadonlyArray1, + PyReadonlyArray2, +}; +use ppvm_lindblad::{codes_from_word, word_from_codes}; +use ppvm_pauli_sum::symmetry as core_sym; +use pyo3::{exceptions::PyValueError, prelude::*}; + +type PyPauliMap<'py> = (Bound<'py, PyArray2>, Bound<'py, PyArray1>); +type PyPauliMapComplex<'py> = (Bound<'py, PyArray2>, Bound<'py, PyArray1>); + +/// A finite abelian symmetry group acting on qubit positions by +/// permutations. Use this to merge translation-equivalent Pauli strings +/// in either the `Lindbladian.pc_step_arr` basis or the `PauliSum` +/// dictionary, reducing per-step memory by up to `|G|×`. +/// +/// Build via the static methods: +/// - `TranslationGroup.chain_1d(n)` — 1D chain of `n` sites with PBC. +/// - `TranslationGroup.torus_2d(lx, ly)` — 2D torus; qubit `(i, j)` at +/// index `j*lx + i`. +/// - `TranslationGroup.torus_3d(lx, ly, lz)` — 3D torus; qubit +/// `(i, j, k)` at index `k*lx*ly + j*lx + i`. +/// - `TranslationGroup.ladder(l, n_legs)` — `n_legs`-leg ladder of `l` +/// sites, translation only along chain direction; qubit `(leg, j)` at +/// index `leg*l + j`. +/// - `TranslationGroup.from_generators(n_qubits, perms, orders)` — +/// arbitrary list of generator permutations + cyclic orders. +#[pyclass(frozen)] +pub struct TranslationGroup { + pub(crate) inner: core_sym::TranslationGroup, +} + +impl TranslationGroup { + /// Accessor for the underlying [`ppvm_pauli_sum::symmetry::TranslationGroup`]. + /// Used by other crate-internal modules (e.g. the PauliSum interface + /// macro) to call into the core merging API. + pub fn core(&self) -> &core_sym::TranslationGroup { + &self.inner + } +} + +#[pymethods] +impl TranslationGroup { + #[staticmethod] + pub fn chain_1d(n: usize) -> Self { + Self { + inner: core_sym::TranslationGroup::chain_1d(n), + } + } + + #[staticmethod] + pub fn torus_2d(lx: usize, ly: usize) -> Self { + Self { + inner: core_sym::TranslationGroup::torus_2d(lx, ly), + } + } + + #[staticmethod] + pub fn torus_3d(lx: usize, ly: usize, lz: usize) -> Self { + Self { + inner: core_sym::TranslationGroup::torus_3d(lx, ly, lz), + } + } + + #[staticmethod] + pub fn ladder(l: usize, n_legs: usize) -> Self { + Self { + inner: core_sym::TranslationGroup::ladder(l, n_legs), + } + } + + #[staticmethod] + pub fn from_generators( + n_qubits: usize, + perms: Vec>, + orders: Vec, + ) -> PyResult { + if perms.len() != orders.len() { + return Err(PyValueError::new_err(format!( + "perms ({} generators) and orders ({}) must have the same length", + perms.len(), + orders.len() + ))); + } + for (g, perm) in perms.iter().enumerate() { + if perm.len() != n_qubits { + return Err(PyValueError::new_err(format!( + "generator {g}: permutation length {} != n_qubits {n_qubits}", + perm.len() + ))); + } + let mut seen = vec![false; n_qubits]; + for &p in perm { + let p = p as usize; + if p >= n_qubits { + return Err(PyValueError::new_err(format!( + "generator {g}: target {p} out of range [0, {n_qubits})" + ))); + } + if seen[p] { + return Err(PyValueError::new_err(format!( + "generator {g}: not a permutation (duplicate target {p})" + ))); + } + seen[p] = true; + } + } + Ok(Self { + inner: core_sym::TranslationGroup::from_generators(n_qubits, perms, orders), + }) + } + + /// Number of qubits this group acts on. + #[getter] + pub fn n_qubits(&self) -> usize { + self.inner.n_qubits() + } + + /// Number of generators (rank as an abelian product group). + #[getter] + pub fn n_generators(&self) -> usize { + self.inner.n_generators() + } + + /// Total group order: product of generator orders. + #[getter] + pub fn order(&self) -> usize { + self.inner.order() + } + + /// Return the canonical (lex-min) orbit representative of `pauli`. + /// `pauli` is a length-`n_qubits` uint8 array with the encoding + /// `0=I, 1=X, 2=Z, 3=Y`. Result is the same shape. + pub fn canonicalize<'py>( + &self, + py: Python<'py>, + pauli: PyReadonlyArray1<'py, u8>, + ) -> PyResult>> { + let codes = pauli.as_slice()?; + if codes.len() != self.inner.n_qubits() { + return Err(PyValueError::new_err(format!( + "pauli has length {} but group expects {} qubits", + codes.len(), + self.inner.n_qubits() + ))); + } + let w = word_from_codes(codes).map_err(|e| PyValueError::new_err(e.to_string()))?; + let canon = self.inner.canonicalize(&w); + let mut out = vec![0u8; codes.len()]; + codes_from_word(&canon, &mut out); + Ok(out.into_pyarray(py)) + } +} + +/// Phase-aware merge of a complex-coefficient `(basis_arr, coeffs)` +/// Pauli sum into orbit-rep form, projected onto momentum sector +/// `momentum`. +/// +/// `momentum` is a length-`group.n_generators` integer array of mode +/// indices; the wavenumber along generator `g` is +/// `2π · momentum[g] / group.generator_order(g)`. Use `momentum=[0, …]` +/// for the trivial (k=0) sector — equivalent to plain merging modulo +/// the 1/|G| normalization the complex merge applies. +/// +/// If the input is **not** in sector `momentum`, the projection +/// silently throws away the other components. Use +/// [`check_momentum_sector_arr`] beforehand to validate. +#[pyfunction] +pub fn canonicalize_basis_arr_complex<'py>( + py: Python<'py>, + basis: PyReadonlyArray2<'py, u8>, + coeffs: PyReadonlyArray1<'py, Complex64>, + group: &TranslationGroup, + momentum: PyReadonlyArray1<'py, i32>, +) -> PyResult> { + let basis_view = basis.as_array(); + let n_q = group.inner.n_qubits(); + if basis_view.shape().get(1).copied() != Some(n_q) { + return Err(PyValueError::new_err(format!( + "basis has {} qubits per row but group acts on {n_q}", + basis_view.shape().get(1).copied().unwrap_or(0) + ))); + } + let n = basis_view.shape()[0]; + let coeffs_slice = coeffs.as_slice()?; + if coeffs_slice.len() != n { + return Err(PyValueError::new_err(format!( + "coeffs has length {} but basis has {} rows", + coeffs_slice.len(), + n + ))); + } + let k_slice = momentum.as_slice()?; + if k_slice.len() != group.inner.n_generators() { + return Err(PyValueError::new_err(format!( + "momentum has {} entries but group has {} generators", + k_slice.len(), + group.inner.n_generators() + ))); + } + let mut basis_words = crate::lindblad::decode_basis(&basis_view, n_q)?; + let mut coeffs_vec: Vec> = coeffs_slice + .iter() + .map(|c| Complex::new(c.re, c.im)) + .collect(); + + core_sym::canonicalize_pauli_sum_complex( + &mut basis_words, + &mut coeffs_vec, + &group.inner, + k_slice, + ); + + let m = basis_words.len(); + let mut out_basis = vec![0u8; m * n_q]; + for (i, w) in basis_words.iter().enumerate() { + codes_from_word(w, &mut out_basis[i * n_q..(i + 1) * n_q]); + } + let out_coeffs: Vec = + coeffs_vec.iter().map(|c| Complex64::new(c.re, c.im)).collect(); + let basis_arr = out_basis + .into_pyarray(py) + .reshape([m, n_q]) + .map_err(|e| PyValueError::new_err(format!("reshape failed: {e}")))?; + Ok((basis_arr, out_coeffs.into_pyarray(py))) +} + +/// Verify that a `(basis_arr, complex_coeffs)` Pauli sum lies in the +/// momentum sector `momentum` under `group`. Returns `None` on pass, +/// raises a `ValueError` with diagnostic info on fail. +/// +/// `tol` is the relative tolerance on coefficient comparison; default +/// `1e-8`. +#[pyfunction] +#[pyo3(signature = (basis, coeffs, group, momentum, tol = 1e-8))] +pub fn check_momentum_sector_arr<'py>( + basis: PyReadonlyArray2<'py, u8>, + coeffs: PyReadonlyArray1<'py, Complex64>, + group: &TranslationGroup, + momentum: PyReadonlyArray1<'py, i32>, + tol: f64, +) -> PyResult<()> { + let basis_view = basis.as_array(); + let n_q = group.inner.n_qubits(); + if basis_view.shape().get(1).copied() != Some(n_q) { + return Err(PyValueError::new_err(format!( + "basis has {} qubits per row but group acts on {n_q}", + basis_view.shape().get(1).copied().unwrap_or(0) + ))); + } + let coeffs_slice = coeffs.as_slice()?; + let k_slice = momentum.as_slice()?; + let basis_words = crate::lindblad::decode_basis(&basis_view, n_q)?; + let coeffs_vec: Vec> = coeffs_slice + .iter() + .map(|c| Complex::new(c.re, c.im)) + .collect(); + core_sym::check_momentum_sector(&basis_words, &coeffs_vec, &group.inner, k_slice, tol) + .map_err(|e| PyValueError::new_err(format!("{e}"))) +} + +/// Merge a `(basis_arr, coeffs)` Pauli sum (the representation used by +/// `Lindbladian.pc_step_arr`) into orbit-representative form. +/// Each row of `basis_arr` is replaced by its canonical +/// representative; coefficients of rows collapsing to the same rep are +/// summed. +/// +/// Returns `(merged_basis_arr, merged_coeffs)`. Output length ≤ input +/// length. +/// +/// For dynamics that commute with `group` and initial states that are +/// `group`-invariant, this preserves all `group`-invariant expectation +/// values (Theorem 1 of Teng et al., arXiv:2512.12094). +#[pyfunction] +pub fn canonicalize_basis_arr<'py>( + py: Python<'py>, + basis: PyReadonlyArray2<'py, u8>, + coeffs: PyReadonlyArray1<'py, f64>, + group: &TranslationGroup, +) -> PyResult> { + let basis_view = basis.as_array(); + let n_q = group.inner.n_qubits(); + if basis_view.shape().get(1).copied() != Some(n_q) { + return Err(PyValueError::new_err(format!( + "basis has {} qubits per row but group acts on {n_q}", + basis_view.shape().get(1).copied().unwrap_or(0) + ))); + } + let n = basis_view.shape()[0]; + let coeffs_slice = coeffs.as_slice()?; + if coeffs_slice.len() != n { + return Err(PyValueError::new_err(format!( + "coeffs has length {} but basis has {} rows", + coeffs_slice.len(), + n + ))); + } + + let mut basis_words = crate::lindblad::decode_basis(&basis_view, n_q)?; + let mut coeffs_vec = coeffs_slice.to_vec(); + + core_sym::canonicalize_pauli_sum(&mut basis_words, &mut coeffs_vec, &group.inner); + + // Re-encode. + let m = basis_words.len(); + let mut out_basis = vec![0u8; m * n_q]; + for (i, w) in basis_words.iter().enumerate() { + codes_from_word(w, &mut out_basis[i * n_q..(i + 1) * n_q]); + } + let basis_arr = out_basis + .into_pyarray(py) + .reshape([m, n_q]) + .map_err(|e| PyValueError::new_err(format!("reshape failed: {e}")))?; + Ok((basis_arr, coeffs_vec.into_pyarray(py))) +} diff --git a/ppvm-python/src/ppvm/lindblad.py b/ppvm-python/src/ppvm/lindblad.py index a31ed4e5f..82dfcd769 100644 --- a/ppvm-python/src/ppvm/lindblad.py +++ b/ppvm-python/src/ppvm/lindblad.py @@ -9,7 +9,8 @@ adaptive Heisenberg-picture evolution: - ``pc_step(...)`` / ``pc_step_arr(...)``: one adaptive predictor-corrector - step ``O ← exp(dt·L*) O`` + step ``O ← exp(dt·L*) O``; ``pc_step_orbit_rep(...)`` is the + translation-symmetric (momentum-sector) variant - ``action(p)`` / ``action_arr(p)``: L*(p) for one Pauli string p - ``leakage(basis, coeffs)`` / ``leakage_arr(...)``: off-basis component of L*(Σ c_j p_j), driving basis expansion @@ -290,6 +291,63 @@ def pc_step_arr( None if tau_add is None else float(tau_add), ) + def pc_step_orbit_rep( + self, + basis_arr: np.ndarray, + coeffs: np.ndarray, + dt: float, + max_basis: int, + group, + momentum: np.ndarray, + drop_tol: float = 1e-12, + protected_arr: np.ndarray | None = None, + canonicalize_first: bool = False, + admit_basis: int | None = None, + tau_add: float | None = None, + ) -> tuple[np.ndarray, np.ndarray]: + """Per-step orbit-representative pc evolution. + + State lives entirely in orbit-rep form throughout: ``basis_arr`` + contains only canonical translation-orbit representatives, + ``coeffs`` are complex, and the action is phase-aware. The basis + is ~``|group|×`` smaller than the equivalent full-basis complex + evolution, and the reduction persists across every step. + + Truncation. ``max_basis`` is a hard rank cap on the live orbit-rep + basis: enrichment adds at most ``max_basis - len(basis)`` of the + largest leakage reps, and the post-step basis is trimmed to the + top-``max_basis`` reps by ``|c|`` (``protected`` reps always kept). + Pass a large value (e.g. ``10_000_000``) for the near-exact, + uncapped case. ``drop_tol`` additionally prunes reps whose absolute + coefficient is below the threshold after the corrector. + + ``admit_basis``, when set (>= ``max_basis``), bounds the enriched + working set instead of ``max_basis``: the step may hold up to + ``admit_basis`` reps transiently and the final truncation keeps the + top-``max_basis`` by evolved ``|c|`` over the whole union — the + displacement scheme, matching the real-space ``pc_step_arr``. + + ``basis_arr`` is assumed to contain canonical reps only. Pass + ``canonicalize_first=True`` to rewrite each row to its canonical + rep on entry (coefficients unchanged). + """ + n = self.n_qubits + if protected_arr is None: + protected_arr = np.zeros((0, n), dtype=np.uint8) + return self._spec.pc_step_orbit_rep( + np.ascontiguousarray(basis_arr, dtype=np.uint8), + np.ascontiguousarray(coeffs, dtype=np.complex128), + float(dt), + int(max_basis), + group, + np.ascontiguousarray(momentum, dtype=np.int32), + float(drop_tol), + np.ascontiguousarray(protected_arr, dtype=np.uint8), + bool(canonicalize_first), + None if admit_basis is None else int(admit_basis), + None if tau_add is None else float(tau_add), + ) + def pc_step( self, basis: Sequence[str], diff --git a/ppvm-python/src/ppvm/paulisum.py b/ppvm-python/src/ppvm/paulisum.py index dca573d8d..b367a7244 100644 --- a/ppvm-python/src/ppvm/paulisum.py +++ b/ppvm-python/src/ppvm/paulisum.py @@ -386,6 +386,50 @@ def trace(self, pattern: str) -> float: """ return self._interface.trace(pattern) + def symmetry_merge(self, group) -> None: + """Merge entries into orbit-representative form under a translation group. + + Each Pauli word in the sum is replaced by its canonical (lex-min) + representative under the action of ``group``; coefficients of words + that collapse to the same representative are summed. Entry count + reduces by up to ``|group|×`` for translation-invariant operators. + + For a translation-invariant dynamics that you apply between + merging steps, this preserves all ``group``-invariant expectation + values (Theorem 1 of Teng et al., arXiv:2512.12094). Plain + real-coefficient merge — handles the trivial (``k=0``) momentum + sector only. + + Args: + group: A `ppvm._core.TranslationGroup` + (use ``TranslationGroup.chain_1d(n)``, ``.torus_2d``, + ``.torus_3d``, ``.ladder``, or ``.from_generators``). + """ + self._interface.symmetry_merge(group) + + def momentum_merge(self, other: "PauliSum", group, momentum) -> None: + """Phase-aware (momentum-sector) merge for a complex operator stored + as a *real pair*: ``self`` is the real part and ``other`` the + imaginary part of ``O = self + i·other``. Both are overwritten in + place with the orbit-representative form projected onto momentum + sector ``momentum``. + + Generalizes `symmetry_merge` to non-trivial momentum sectors + (``k != 0``) while keeping real coefficients on both PauliSums — the + only complex arithmetic is the internal character-weighted fold. + ``self`` and ``other`` must be distinct objects with the same qubit + count. Exact after a translation-covariant gate layer; under a + generic Trotter step it carries the same ``O(dt^{p+1})`` equivariance + error as the ``k=0`` merge. + + Args: + other: the PauliSum holding the imaginary component (modified in place). + group: a `ppvm._core.TranslationGroup`. + momentum: sequence of integer modes, one per group generator + (e.g. ``[k]`` for a 1D chain; ``[0, ...]`` is the trivial sector). + """ + self._interface.momentum_merge(other._interface, group, list(momentum)) + def amplitude_damping(self, addr0: int, gamma: float, *, truncate: bool = True): """Apply an amplitude-damping channel. diff --git a/ppvm-python/test/test_momentum_merge.py b/ppvm-python/test/test_momentum_merge.py new file mode 100644 index 000000000..de3302ffa --- /dev/null +++ b/ppvm-python/test/test_momentum_merge.py @@ -0,0 +1,199 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 + +"""Tests for momentum-sector (k != 0) symmetry merging of real PauliSum pairs. + +A complex operator O = O_re + i·O_im is carried as a pair of real PauliSums. +``PauliSum.momentum_merge`` folds the pair onto translation-orbit +representatives in momentum sector k, generalizing ``symmetry_merge`` (k=0). + +These checks compare against *exact* references — the projector definition, +idempotency, and exact diagonalization of the dynamics — NOT against any +other propagation scheme. +""" +import cmath +import math + +import numpy as np +import pytest + +from ppvm import PauliSum +from ppvm._core import TranslationGroup + +# ── dense Pauli helpers (exact references) ─────────────────────────────────── +_I = np.eye(2, dtype=complex) +_X = np.array([[0, 1], [1, 0]], dtype=complex) +_Y = np.array([[0, -1j], [1j, 0]], dtype=complex) +_Z = np.array([[1, 0], [0, -1]], dtype=complex) +_P = {"I": _I, "X": _X, "Y": _Y, "Z": _Z} + + +def dense(pauli_str): + m = np.array([[1]], dtype=complex) + for ch in pauli_str: + m = np.kron(m, _P[ch]) + return m + + +def zstr(n, q): + return "".join("Z" if i == q else "I" for i in range(n)) + + +def chain_bonds(n): + return [(i, (i + 1) % n, 1.0) for i in range(n)] + + +# ── helpers shared with the k-resolved Trotter driver ──────────────────────── +def _seed_pair(n, k): + a = np.arange(n) + re = np.cos(2 * np.pi * k * a / n) + im = -np.sin(2 * np.pi * k * a / n) # e^{-2πi k a/n} = cos - i sin + Z = [zstr(n, q) for q in range(n)] + PA = PauliSum.new(n, [(Z[q], float(re[q])) for q in range(n)], + min_abs_coeff=0.0, max_pauli_weight=n) + PB = PauliSum.new(n, [(Z[q], float(im[q])) for q in range(n)], + min_abs_coeff=0.0, max_pauli_weight=n) + return PA, PB + + +def _to_complex_dict(PA, PB): + d = {} + for s, c in PA.terms: + d[s] = d.get(s, 0j) + c + for s, c in PB.terms: + d[s] = d.get(s, 0j) + 1j * c + return {s: v for s, v in d.items() if v != 0j} + + +def _ovl(sA, sB, oA, oB): + re = sA.overlap(oA) + sB.overlap(oB) + im = sA.overlap(oB) - sB.overlap(oA) + return complex(re, im) + + +# ============================================================================= +# 1. The merge is an exact sector projector: idempotent, and it leaves a +# genuine momentum-k eigenoperator unchanged. +# ============================================================================= +@pytest.mark.parametrize("k", [0, 1, 2, 3]) +def test_momentum_merge_idempotent(k): + n = 4 + g = TranslationGroup.chain_1d(n) + PA, PB = _seed_pair(n, k) # S^z_k is exactly in sector k + PA.momentum_merge(PB, g, [k]) + once = _to_complex_dict(PA, PB) + PA.momentum_merge(PB, g, [k]) # merging again must be a no-op + twice = _to_complex_dict(PA, PB) + keys = set(once) | set(twice) + assert max(abs(once.get(x, 0j) - twice.get(x, 0j)) for x in keys) < 1e-12 + + +def test_momentum_merge_projects_out_other_sectors(): + """Merging a pure sector-k operator in sector k' != k gives ~zero.""" + n = 4 + g = TranslationGroup.chain_1d(n) + PA, PB = _seed_pair(n, 1) # operator lives in k=1 + PA.momentum_merge(PB, g, [2]) # project onto k=2 + d = _to_complex_dict(PA, PB) + assert all(abs(v) < 1e-12 for v in d.values()), d + + +# ============================================================================= +# 2. End-to-end: k-resolved, symmetry-compressed Trotter reproduces the +# EXACT (dense-diagonalization) operator autocorrelator as dt -> 0. +# ============================================================================= +def _ed_autocorr(n, bonds, k, ts): + """C_k(t) = Tr[O0^dagger O(t)] / Tr[O0^dagger O0], O0 = S^z_k, exact.""" + H = np.zeros((2 ** n, 2 ** n), dtype=complex) + for (i, j, J) in bonds: + for q in "XY": + s = ["I"] * n + s[i] = q + s[j] = q + H += J * dense("".join(s)) + O0 = np.zeros((2 ** n, 2 ** n), dtype=complex) + for a in range(n): + O0 += cmath.exp(-2j * math.pi * k * a / n) * dense(zstr(n, a)) + E, V = np.linalg.eigh(H) + out = [] + with np.errstate(all="ignore"): # silence spurious macOS-Accelerate matmul warnings + norm = np.trace(O0.conj().T @ O0).real + for t in ts: + U = (V * np.exp(-1j * E * t)) @ V.conj().T + Ot = U.conj().T @ O0 @ U + out.append(np.trace(O0.conj().T @ Ot) / norm) + return np.array(out) + + +def _ctrotter_autocorr(n, bonds, k, dt, steps): + g = TranslationGroup.chain_1d(n) + PA, PB = _seed_pair(n, k) + PA.momentum_merge(PB, g, [k]) + refA, refB = PA.copy(), PB.copy() + C0 = _ovl(refA, refB, PA, PB) + out = [1.0 + 0j] + for _ in range(steps): + for (i, j, J) in bonds: # Strang: forward then reversed + PA.rxx(i, j, J * dt, truncate=False); PA.ryy(i, j, J * dt, truncate=False) + PB.rxx(i, j, J * dt, truncate=False); PB.ryy(i, j, J * dt, truncate=False) + for (i, j, J) in reversed(bonds): + PA.rxx(i, j, J * dt, truncate=False); PA.ryy(i, j, J * dt, truncate=False) + PB.rxx(i, j, J * dt, truncate=False); PB.ryy(i, j, J * dt, truncate=False) + PA.momentum_merge(PB, g, [k]) + out.append(_ovl(refA, refB, PA, PB) / C0) + return np.array(out) + + +@pytest.mark.parametrize("k", [0, 1, 2, 3]) +def test_k_resolved_trotter_converges_to_exact(k): + n, T = 4, 0.3 + bonds = chain_bonds(n) + # exact reference at the matching times for two step sizes + err = {} + for dt in (0.04, 0.02): + steps = round(T / dt) + ts = np.arange(steps + 1) * dt + c = _ctrotter_autocorr(n, bonds, k, dt, steps) + ed = _ed_autocorr(n, bonds, k, ts) + err[dt] = np.max(np.abs(c - ed)) + + assert abs(_ctrotter_autocorr(n, bonds, k, 0.02, 1)[0] - 1.0) < 1e-12 # C_k(0)=1 + if k == 0: + # total Z is conserved -> exact in every sector-0 step + assert err[0.02] < 1e-10 + else: + assert err[0.02] < 5e-3 # close to exact at dt=0.02 + assert err[0.02] < err[0.04] # converges toward exact as dt->0 + + +def test_compressed_matches_uncompressed_evolution(): + """Merging must not change observables beyond the O(dt^2) equivariance + error: compressed (merge each step) vs the same gates with no merge.""" + n, k, dt, steps = 4, 2, 0.02, 10 + bonds = chain_bonds(n) + g = TranslationGroup.chain_1d(n) + + # uncompressed: evolve the full real pair, project only at readout + PA, PB = _seed_pair(n, k) + rA, rB = _seed_pair(n, k) + rA.momentum_merge(rB, g, [k]) + C0 = _ovl(rA, rB, *_merged_copy(PA, PB, g, k)) + comp = _ctrotter_autocorr(n, bonds, k, dt, steps) + unc = [] + for _ in range(steps): + for (i, j, J) in bonds: + PA.rxx(i, j, J * dt, truncate=False); PA.ryy(i, j, J * dt, truncate=False) + PB.rxx(i, j, J * dt, truncate=False); PB.ryy(i, j, J * dt, truncate=False) + for (i, j, J) in reversed(bonds): + PA.rxx(i, j, J * dt, truncate=False); PA.ryy(i, j, J * dt, truncate=False) + PB.rxx(i, j, J * dt, truncate=False); PB.ryy(i, j, J * dt, truncate=False) + mA, mB = _merged_copy(PA, PB, g, k) + unc.append(_ovl(rA, rB, mA, mB) / C0) + unc = np.array([1.0 + 0j] + unc) + assert np.max(np.abs(comp - unc)) < 5e-3 # only O(dt^2) equivariance + + +def _merged_copy(PA, PB, g, k): + a, b = PA.copy(), PB.copy() + a.momentum_merge(b, g, [k]) + return a, b From ece442e98e7aad9deb2a9d002e0bfd94d104217f Mon Sep 17 00:00:00 2001 From: alexschuckert Date: Thu, 16 Jul 2026 13:35:12 +0100 Subject: [PATCH 02/13] style: cargo fmt + ruff format Co-Authored-By: Claude Fable 5 --- crates/ppvm-lindblad/src/orbit_rep.rs | 28 ++++++++-- crates/ppvm-python-native/src/interface.rs | 5 +- crates/ppvm-python-native/src/symmetry.rs | 9 ++-- ppvm-python/test/test_momentum_merge.py | 63 +++++++++++++--------- 4 files changed, 66 insertions(+), 39 deletions(-) diff --git a/crates/ppvm-lindblad/src/orbit_rep.rs b/crates/ppvm-lindblad/src/orbit_rep.rs index 7bf7a27c6..c72e3da48 100644 --- a/crates/ppvm-lindblad/src/orbit_rep.rs +++ b/crates/ppvm-lindblad/src/orbit_rep.rs @@ -159,7 +159,10 @@ pub(crate) fn expm_apply_orbit_rep_cached( let (m_star, s) = crate::expm::select_ms(dt.abs() * onenorm); let mut v = coeffs.to_vec(); - let op = CscOp { cols: &cols, dim: n }; + let op = CscOp { + cols: &cols, + dim: n, + }; let e = ExpmOp::from_parts( op, Complex::new(dt, 0.0), @@ -310,7 +313,13 @@ pub fn pc_step_orbit_rep( k_modes: &[i32], cfg: &crate::PcStepConfig, ) -> Result<(), Error> { - let crate::PcStepConfig { max_basis, admit_basis, drop_tol, tau_add, .. } = *cfg; + let crate::PcStepConfig { + max_basis, + admit_basis, + drop_tol, + tau_add, + .. + } = *cfg; // Admission bound, mirroring the real-space `pc_step`: enrichment may // grow the live basis to `admit` >= `max_basis`; the final // `cap_basis_complex` keeps the top-`max_basis` reps by evolved |coeff| @@ -329,8 +338,15 @@ pub fn pc_step_orbit_rep( // once via `build_orbit_rep_cols`). let coeffs_predict = expm_apply_orbit_rep_cached(spec, basis, group, k_modes, dt, coeffs); // 3. Second-hop leakage from predicted state. - let mut leak2 = - leakage_orbit_rep(spec, basis, &coeffs_predict, protected, group, k_modes, admit)?; + let mut leak2 = leakage_orbit_rep( + spec, + basis, + &coeffs_predict, + protected, + group, + k_modes, + admit, + )?; drop(coeffs_predict); if tau_add > 0.0 { leak2.retain(|(_, c)| c.norm() > tau_add); @@ -401,7 +417,9 @@ fn cap_basis_complex( return; } else { let k = slots - 1; - mags.select_nth_unstable_by(k, |a, b| b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal)); + mags.select_nth_unstable_by(k, |a, b| { + b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal) + }); mags[k] }; let mut write = 0; diff --git a/crates/ppvm-python-native/src/interface.rs b/crates/ppvm-python-native/src/interface.rs index ec690ae80..54813f130 100644 --- a/crates/ppvm-python-native/src/interface.rs +++ b/crates/ppvm-python-native/src/interface.rs @@ -78,10 +78,7 @@ macro_rules! create_interface_symmetry_methods { group.core().n_qubits(), ))); } - ppvm_pauli_sum::symmetry::symmetry_merge_pauli_sum( - &mut self.inner, - group.core(), - ); + ppvm_pauli_sum::symmetry::symmetry_merge_pauli_sum(&mut self.inner, group.core()); Ok(()) } diff --git a/crates/ppvm-python-native/src/symmetry.rs b/crates/ppvm-python-native/src/symmetry.rs index d3c69b6e9..880851468 100644 --- a/crates/ppvm-python-native/src/symmetry.rs +++ b/crates/ppvm-python-native/src/symmetry.rs @@ -12,8 +12,7 @@ use num::Complex; use numpy::{ - Complex64, IntoPyArray, PyArray1, PyArray2, PyArrayMethods, PyReadonlyArray1, - PyReadonlyArray2, + Complex64, IntoPyArray, PyArray1, PyArray2, PyArrayMethods, PyReadonlyArray1, PyReadonlyArray2, }; use ppvm_lindblad::{codes_from_word, word_from_codes}; use ppvm_pauli_sum::symmetry as core_sym; @@ -229,8 +228,10 @@ pub fn canonicalize_basis_arr_complex<'py>( for (i, w) in basis_words.iter().enumerate() { codes_from_word(w, &mut out_basis[i * n_q..(i + 1) * n_q]); } - let out_coeffs: Vec = - coeffs_vec.iter().map(|c| Complex64::new(c.re, c.im)).collect(); + let out_coeffs: Vec = coeffs_vec + .iter() + .map(|c| Complex64::new(c.re, c.im)) + .collect(); let basis_arr = out_basis .into_pyarray(py) .reshape([m, n_q]) diff --git a/ppvm-python/test/test_momentum_merge.py b/ppvm-python/test/test_momentum_merge.py index de3302ffa..8566fe2f1 100644 --- a/ppvm-python/test/test_momentum_merge.py +++ b/ppvm-python/test/test_momentum_merge.py @@ -11,6 +11,7 @@ idempotency, and exact diagonalization of the dynamics — NOT against any other propagation scheme. """ + import cmath import math @@ -49,10 +50,12 @@ def _seed_pair(n, k): re = np.cos(2 * np.pi * k * a / n) im = -np.sin(2 * np.pi * k * a / n) # e^{-2πi k a/n} = cos - i sin Z = [zstr(n, q) for q in range(n)] - PA = PauliSum.new(n, [(Z[q], float(re[q])) for q in range(n)], - min_abs_coeff=0.0, max_pauli_weight=n) - PB = PauliSum.new(n, [(Z[q], float(im[q])) for q in range(n)], - min_abs_coeff=0.0, max_pauli_weight=n) + PA = PauliSum.new( + n, [(Z[q], float(re[q])) for q in range(n)], min_abs_coeff=0.0, max_pauli_weight=n + ) + PB = PauliSum.new( + n, [(Z[q], float(im[q])) for q in range(n)], min_abs_coeff=0.0, max_pauli_weight=n + ) return PA, PB @@ -79,10 +82,10 @@ def _ovl(sA, sB, oA, oB): def test_momentum_merge_idempotent(k): n = 4 g = TranslationGroup.chain_1d(n) - PA, PB = _seed_pair(n, k) # S^z_k is exactly in sector k + PA, PB = _seed_pair(n, k) # S^z_k is exactly in sector k PA.momentum_merge(PB, g, [k]) once = _to_complex_dict(PA, PB) - PA.momentum_merge(PB, g, [k]) # merging again must be a no-op + PA.momentum_merge(PB, g, [k]) # merging again must be a no-op twice = _to_complex_dict(PA, PB) keys = set(once) | set(twice) assert max(abs(once.get(x, 0j) - twice.get(x, 0j)) for x in keys) < 1e-12 @@ -92,8 +95,8 @@ def test_momentum_merge_projects_out_other_sectors(): """Merging a pure sector-k operator in sector k' != k gives ~zero.""" n = 4 g = TranslationGroup.chain_1d(n) - PA, PB = _seed_pair(n, 1) # operator lives in k=1 - PA.momentum_merge(PB, g, [2]) # project onto k=2 + PA, PB = _seed_pair(n, 1) # operator lives in k=1 + PA.momentum_merge(PB, g, [2]) # project onto k=2 d = _to_complex_dict(PA, PB) assert all(abs(v) < 1e-12 for v in d.values()), d @@ -104,14 +107,14 @@ def test_momentum_merge_projects_out_other_sectors(): # ============================================================================= def _ed_autocorr(n, bonds, k, ts): """C_k(t) = Tr[O0^dagger O(t)] / Tr[O0^dagger O0], O0 = S^z_k, exact.""" - H = np.zeros((2 ** n, 2 ** n), dtype=complex) - for (i, j, J) in bonds: + H = np.zeros((2**n, 2**n), dtype=complex) + for i, j, J in bonds: for q in "XY": s = ["I"] * n s[i] = q s[j] = q H += J * dense("".join(s)) - O0 = np.zeros((2 ** n, 2 ** n), dtype=complex) + O0 = np.zeros((2**n, 2**n), dtype=complex) for a in range(n): O0 += cmath.exp(-2j * math.pi * k * a / n) * dense(zstr(n, a)) E, V = np.linalg.eigh(H) @@ -133,12 +136,16 @@ def _ctrotter_autocorr(n, bonds, k, dt, steps): C0 = _ovl(refA, refB, PA, PB) out = [1.0 + 0j] for _ in range(steps): - for (i, j, J) in bonds: # Strang: forward then reversed - PA.rxx(i, j, J * dt, truncate=False); PA.ryy(i, j, J * dt, truncate=False) - PB.rxx(i, j, J * dt, truncate=False); PB.ryy(i, j, J * dt, truncate=False) - for (i, j, J) in reversed(bonds): - PA.rxx(i, j, J * dt, truncate=False); PA.ryy(i, j, J * dt, truncate=False) - PB.rxx(i, j, J * dt, truncate=False); PB.ryy(i, j, J * dt, truncate=False) + for i, j, J in bonds: # Strang: forward then reversed + PA.rxx(i, j, J * dt, truncate=False) + PA.ryy(i, j, J * dt, truncate=False) + PB.rxx(i, j, J * dt, truncate=False) + PB.ryy(i, j, J * dt, truncate=False) + for i, j, J in reversed(bonds): + PA.rxx(i, j, J * dt, truncate=False) + PA.ryy(i, j, J * dt, truncate=False) + PB.rxx(i, j, J * dt, truncate=False) + PB.ryy(i, j, J * dt, truncate=False) PA.momentum_merge(PB, g, [k]) out.append(_ovl(refA, refB, PA, PB) / C0) return np.array(out) @@ -162,8 +169,8 @@ def test_k_resolved_trotter_converges_to_exact(k): # total Z is conserved -> exact in every sector-0 step assert err[0.02] < 1e-10 else: - assert err[0.02] < 5e-3 # close to exact at dt=0.02 - assert err[0.02] < err[0.04] # converges toward exact as dt->0 + assert err[0.02] < 5e-3 # close to exact at dt=0.02 + assert err[0.02] < err[0.04] # converges toward exact as dt->0 def test_compressed_matches_uncompressed_evolution(): @@ -181,16 +188,20 @@ def test_compressed_matches_uncompressed_evolution(): comp = _ctrotter_autocorr(n, bonds, k, dt, steps) unc = [] for _ in range(steps): - for (i, j, J) in bonds: - PA.rxx(i, j, J * dt, truncate=False); PA.ryy(i, j, J * dt, truncate=False) - PB.rxx(i, j, J * dt, truncate=False); PB.ryy(i, j, J * dt, truncate=False) - for (i, j, J) in reversed(bonds): - PA.rxx(i, j, J * dt, truncate=False); PA.ryy(i, j, J * dt, truncate=False) - PB.rxx(i, j, J * dt, truncate=False); PB.ryy(i, j, J * dt, truncate=False) + for i, j, J in bonds: + PA.rxx(i, j, J * dt, truncate=False) + PA.ryy(i, j, J * dt, truncate=False) + PB.rxx(i, j, J * dt, truncate=False) + PB.ryy(i, j, J * dt, truncate=False) + for i, j, J in reversed(bonds): + PA.rxx(i, j, J * dt, truncate=False) + PA.ryy(i, j, J * dt, truncate=False) + PB.rxx(i, j, J * dt, truncate=False) + PB.ryy(i, j, J * dt, truncate=False) mA, mB = _merged_copy(PA, PB, g, k) unc.append(_ovl(rA, rB, mA, mB) / C0) unc = np.array([1.0 + 0j] + unc) - assert np.max(np.abs(comp - unc)) < 5e-3 # only O(dt^2) equivariance + assert np.max(np.abs(comp - unc)) < 5e-3 # only O(dt^2) equivariance def _merged_copy(PA, PB, g, k): From d769a94dd77869ddacc5ff158538d9e426a790ac Mon Sep 17 00:00:00 2001 From: alexschuckert Date: Thu, 16 Jul 2026 13:52:36 +0100 Subject: [PATCH 03/13] style: ruff RUF005 (unpack instead of concatenation) Co-Authored-By: Claude Fable 5 --- ppvm-python/test/test_momentum_merge.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ppvm-python/test/test_momentum_merge.py b/ppvm-python/test/test_momentum_merge.py index 8566fe2f1..0546deb11 100644 --- a/ppvm-python/test/test_momentum_merge.py +++ b/ppvm-python/test/test_momentum_merge.py @@ -200,7 +200,7 @@ def test_compressed_matches_uncompressed_evolution(): PB.ryy(i, j, J * dt, truncate=False) mA, mB = _merged_copy(PA, PB, g, k) unc.append(_ovl(rA, rB, mA, mB) / C0) - unc = np.array([1.0 + 0j] + unc) + unc = np.array([1.0 + 0j, *unc]) assert np.max(np.abs(comp - unc)) < 5e-3 # only O(dt^2) equivariance From 4a58ab3ce5006191670ac171ffec5e74f731c71d Mon Sep 17 00:00:00 2001 From: alexschuckert Date: Thu, 16 Jul 2026 13:57:00 +0100 Subject: [PATCH 04/13] types: TranslationGroup / merge / orbit-step stubs in _core.pyi Co-Authored-By: Claude Fable 5 --- ppvm-python/src/ppvm/_core.pyi | 60 ++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/ppvm-python/src/ppvm/_core.pyi b/ppvm-python/src/ppvm/_core.pyi index 4b2cfac7f..b0350fc55 100644 --- a/ppvm-python/src/ppvm/_core.pyi +++ b/ppvm-python/src/ppvm/_core.pyi @@ -62,6 +62,14 @@ class _PauliSumBase: def terms(self) -> list[tuple[str, float]]: ... def weights(self) -> list[tuple[str, int]]: ... def current_max_weight(self) -> int: ... + # Only on non-loss variants (see create_interface_symmetry_methods). + def symmetry_merge(self, group: TranslationGroup) -> None: ... + def momentum_merge( + self, + other: _PauliSumBase, + group: TranslationGroup, + momentum: list[int], + ) -> None: ... class _PauliSumLossBase(_PauliSumBase): def loss_channel(self, addr0: int, p: float, truncate: bool = True) -> None: ... @@ -392,4 +400,56 @@ class LindbladSpec: admit_basis: int | None = None, tau_add: float | None = None, ) -> tuple[tuple[np.ndarray, np.ndarray], dict[str, int]]: ... + def pc_step_orbit_rep( + self, + basis: np.ndarray, + coeffs: np.ndarray, + dt: float, + max_basis: int, + group: TranslationGroup, + momentum: np.ndarray, + drop_tol: float = 0.0, + protected: np.ndarray | None = None, + canonicalize_first: bool = False, + admit_basis: int | None = None, + tau_add: float | None = None, + ) -> tuple[np.ndarray, np.ndarray]: ... def generator(self, basis: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]: ... + +class TranslationGroup: + @staticmethod + def chain_1d(n: int) -> TranslationGroup: ... + @staticmethod + def torus_2d(lx: int, ly: int) -> TranslationGroup: ... + @staticmethod + def torus_3d(lx: int, ly: int, lz: int) -> TranslationGroup: ... + @staticmethod + def ladder(l: int, n_legs: int) -> TranslationGroup: ... + @staticmethod + def from_generators( + n_qubits: int, perms: list[list[int]], orders: list[int] + ) -> TranslationGroup: ... + @property + def n_qubits(self) -> int: ... + @property + def n_generators(self) -> int: ... + @property + def order(self) -> int: ... + def canonicalize(self, pauli: np.ndarray) -> np.ndarray: ... + +def canonicalize_basis_arr( + basis: np.ndarray, coeffs: np.ndarray, group: TranslationGroup +) -> tuple[np.ndarray, np.ndarray]: ... +def canonicalize_basis_arr_complex( + basis: np.ndarray, + coeffs: np.ndarray, + group: TranslationGroup, + momentum: np.ndarray, +) -> tuple[np.ndarray, np.ndarray]: ... +def check_momentum_sector_arr( + basis: np.ndarray, + coeffs: np.ndarray, + group: TranslationGroup, + momentum: np.ndarray, + tol: float = 1e-8, +) -> None: ... From 070bc753b0df30f6290b2d62415013d7cbc13190 Mon Sep 17 00:00:00 2001 From: alexschuckert Date: Thu, 16 Jul 2026 16:13:01 +0100 Subject: [PATCH 05/13] style: type alias for the orbit-step return (clippy type_complexity) Co-Authored-By: Claude Fable 5 --- crates/ppvm-python-native/src/lindblad.rs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/crates/ppvm-python-native/src/lindblad.rs b/crates/ppvm-python-native/src/lindblad.rs index 07d0e739f..7c7206201 100644 --- a/crates/ppvm-python-native/src/lindblad.rs +++ b/crates/ppvm-python-native/src/lindblad.rs @@ -20,6 +20,7 @@ use ppvm_lindblad::{JumpInput, LindbladSpec as CoreSpec, Word, codes_from_word, use pyo3::{exceptions::PyValueError, prelude::*}; type PyPauliMap<'py> = (Bound<'py, PyArray2>, Bound<'py, PyArray1>); +type PyPauliMapComplex<'py> = (Bound<'py, PyArray2>, Bound<'py, PyArray1>); type PyCoo<'py> = ( Bound<'py, PyArray1>, Bound<'py, PyArray1>, @@ -388,7 +389,7 @@ impl LindbladSpec { canonicalize_first: bool, admit_basis: Option, tau_add: Option, - ) -> PyResult<(Bound<'py, PyArray2>, Bound<'py, PyArray1>)> { + ) -> PyResult> { use num::Complex; use ppvm_lindblad::orbit_rep; From b0d1752e95dad3449f8962d83baec9d1e5a9b92a Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 2 Sep 2026 09:40:40 +0200 Subject: [PATCH 06/13] Fix pre-commit failure --- crates/ppvm-lindblad/src/lib.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/crates/ppvm-lindblad/src/lib.rs b/crates/ppvm-lindblad/src/lib.rs index e99bbb5d0..3af2397dd 100644 --- a/crates/ppvm-lindblad/src/lib.rs +++ b/crates/ppvm-lindblad/src/lib.rs @@ -58,4 +58,3 @@ pub use word::{MAX_QUBITS, Word, codes_from_word, parse_pauli_string, word_from_ #[cfg(test)] mod tests; - From 8d533275ec4eb275dfdb22746f955c6d5c31294e Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 2 Sep 2026 11:34:12 +0200 Subject: [PATCH 07/13] refactor(ppvm-lindblad): dissolve orbit_rep.rs into the spec/basis/step layout MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `orbit_rep.rs` (463 lines) re-implemented one function from each of the crate's existing modules, so every function in it had a near-twin — three of the truncation helpers were verbatim copies with `.abs()` -> `.norm()`. Each function now sits beside its twin, on shared generic code: - new `sector.rs`: `Sector { group, k_modes }` bundles the two values that every phase-aware routine needs together, with one `canonicalize_phase` method replacing the canonicalize/character/multiply idiom at its three call sites. Carries the orbit-rep narrative docs. - new `truncate.rs`: one `prune_basis` / `cap_basis` / `add_leakage_capped` (plus `cap_map_to_room`, `order_by_desc_mag`) generic over the new `scalar::Coeff`, replacing five real/complex duplicate pairs. - `mf_expm.rs`: `build_orbit_rep_cols` + `expm_apply_orbit_rep` beside the real pair, both on a shared `expm_apply_cached` tail (mu, 1-norm, `from_parts`, `apply`). The paths now differ only in the `(m*, s, tol)` selection closure. - `basis.rs`: `leakage_orbit_rep` beside `leakage_complex`. - `step.rs`: `pc_step_orbit_rep` beside `pc_step`. The `Sector` bundle drops `leakage_orbit_rep` 7->6 args and `pc_step_orbit_rep` 8->7, so both `#[allow(clippy::too_many_arguments)]` suppressions are gone; the crate now has zero clippy allows. No behavior change: 348 insertions, 659 deletions. `cargo clippy --workspace --all-targets` is clean, `cargo test --workspace` and the 239 Python tests pass, including the orbit-rep vs full-basis exactness test. Co-Authored-By: Claude Opus 5 (1M context) --- crates/ppvm-lindblad/src/basis.rs | 117 ++++-- crates/ppvm-lindblad/src/lib.rs | 8 +- crates/ppvm-lindblad/src/mf_expm.rs | 206 +++++++--- crates/ppvm-lindblad/src/orbit_rep.rs | 463 ---------------------- crates/ppvm-lindblad/src/scalar.rs | 43 ++ crates/ppvm-lindblad/src/sector.rs | 75 ++++ crates/ppvm-lindblad/src/step.rs | 169 ++++---- crates/ppvm-lindblad/src/tests.rs | 7 +- crates/ppvm-lindblad/src/truncate.rs | 153 +++++++ crates/ppvm-python-native/src/lindblad.rs | 37 +- 10 files changed, 619 insertions(+), 659 deletions(-) delete mode 100644 crates/ppvm-lindblad/src/orbit_rep.rs create mode 100644 crates/ppvm-lindblad/src/scalar.rs create mode 100644 crates/ppvm-lindblad/src/sector.rs create mode 100644 crates/ppvm-lindblad/src/truncate.rs diff --git a/crates/ppvm-lindblad/src/basis.rs b/crates/ppvm-lindblad/src/basis.rs index 9f352d972..941c67b6c 100644 --- a/crates/ppvm-lindblad/src/basis.rs +++ b/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 { @@ -71,17 +77,7 @@ impl LindbladSpec { let protected_set: FxHashMap = 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 = (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 = 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 = 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 = protected.iter().map(|w| (word_hash(w), ())).collect(); - const CHUNK_SIZE: usize = 4096; let n_qubits = self.n_qubits(); let mut merged: FxHashMap> = FxHashMap::default(); for chunk_start in (0..basis.len()).step_by(CHUNK_SIZE) { @@ -253,4 +234,84 @@ 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)` via [`Sector::canonicalize_phase`]. + /// 2. If `r_q` NOT in `basis` and NOT in `protected`: + /// `merged[r_q] += χ_k · v_q · c_r`. + /// + /// 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], + protected: &[Word], + sector: Sector<'_>, + max_basis: usize, + ) -> Result)>, 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> = FxHashMap::default(); + for chunk_indices in order.chunks(CHUNK_SIZE) { + let local: Vec)>> = chunk_indices + .par_iter() + .map_init( + || { + ( + Vec::::with_capacity(n_qubits), + Vec::::with_capacity(128), + FxHashMap::>::with_capacity_and_hasher( + 128, + FxBuildHasher::default(), + ), + ) + }, + |(s1, s2, lm), &i| { + let r = &basis[i]; + let c_r = coeffs[i]; + 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 (r_q, phase) = sector.canonicalize_phase(q); + if !in_basis.contains(&r_q) && !protected_set.contains(&r_q) { + out.push((r_q, phase * *v * c_r)); + } + } + 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()) + } } diff --git a/crates/ppvm-lindblad/src/lib.rs b/crates/ppvm-lindblad/src/lib.rs index 3af2397dd..57dff8aa2 100644 --- a/crates/ppvm-lindblad/src/lib.rs +++ b/crates/ppvm-lindblad/src/lib.rs @@ -38,20 +38,20 @@ 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. pub(crate) mod mf_expm; -/// Per-step orbit-rep evolution under translation symmetry, with a -/// phase-aware complex action. See module docs. -pub mod orbit_rep; - 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}; diff --git a/crates/ppvm-lindblad/src/mf_expm.rs b/crates/ppvm-lindblad/src/mf_expm.rs index f2e5b73f3..491f686d5 100644 --- a/crates/ppvm-lindblad/src/mf_expm.rs +++ b/crates/ppvm-lindblad/src/mf_expm.rs @@ -1,12 +1,14 @@ // SPDX-FileCopyrightText: 2026 The PPVM Authors // SPDX-License-Identifier: Apache-2.0 -//! Matrix-free `exp(dt · L*) · b` for the real (`f64`) path, driven by the -//! external `quspin-expm` crate. +//! Matrix-free `exp(dt · L*) · b`, driven by the external `quspin-expm` +//! crate — for both the real (`f64`) adaptive path and the complex, +//! phase-aware orbit-rep path. //! //! Instead of materialising the in-basis-restricted generator as a CSR, the //! per-column generator action is computed ONCE per expm call (via -//! [`build_mf_cols`]) and reused, CSC-style, across every Krylov/Taylor matvec +//! [`build_mf_cols`] / [`build_orbit_rep_cols`]) and reused, CSC-style, +//! across every Krylov/Taylor matvec //! by [`CscOp`] (a [`quspin_types::LinearOperator`]) fed to //! [`quspin_expm::ExpmOp::from_parts`]. Each matvec is then a cheap CSC //! scatter; the Pauli-commutator action is never recomputed per matvec. @@ -16,17 +18,28 @@ //! `dot_transpose` are never invoked on the single-vector `apply` path; only //! [`LinearOperator::dot`] runs. //! -//! `μ`, the trace, and the exact column 1-norm of `A − μ·I` are computed in -//! the same single action pass as the cache. The `(m, s)` Taylor partition is +//! `μ`, the trace, and the column 1-norm of `A − μ·I` are computed in the +//! same single action pass as the cache, and turned into an `apply` by the +//! shared [`expm_apply_cached`] tail. The `(m, s)` Taylor partition is //! picked with the tolerance-matched tables in [`crate::expm`]: a relaxed //! `tol=1e-6` table when the PC prunes coarsely (`drop_tol ≥ 1e-4`), else the //! double-precision table (keeping the exact-reference test paths bit-exact). +use crate::scalar::Coeff; +use crate::sector::Sector; use crate::{LindbladSpec, Word, build_basis_index, expm}; use fxhash::{FxBuildHasher, FxHashMap}; use num::Complex; use quspin_types::{ExpmComputation, LinearOperator, QuSpinError}; use rayon::prelude::*; +use std::iter::Sum; +use std::ops::{AddAssign, Div, Mul, Sub}; + +/// CSC columns of a cached in-basis action: `cols[c]` = `(row, coeff)`. +type Cols = Vec>; +/// Per-column `(raw, diag)` for the `μ`/1-norm selection: `raw` bounds +/// `Σ_r |M[r,c]|` from above and `diag = M[c,c]`. +type PerCol = Vec<(f64, T)>; /// Per-column in-basis action of the real generator `M`, plus the data the /// `(m, s)`/`μ` selection needs — all from ONE action pass over the basis. @@ -37,16 +50,11 @@ use rayon::prelude::*; /// outputs (in- and out-of-basis, an upper bound on the column 1-norm) and /// `diag` the coefficient of the output Word equal to the input Word. The /// cache is reused by [`CscOp`] across every Krylov/Taylor matvec. -/// CSC columns of the cached in-basis action: `cols[c]` = `(row, coeff)`. -type MfCols = Vec>; -/// Per-column `(raw, diag)` for the `μ`/1-norm selection. -type MfPerCol = Vec<(f64, f64)>; - fn build_mf_cols( spec: &LindbladSpec, basis: &[Word], index: &FxHashMap, -) -> (MfCols, MfPerCol) { +) -> (Cols, PerCol) { basis .par_iter() .map_init( @@ -80,6 +88,63 @@ fn build_mf_cols( .unzip() } +/// Per-column **phase-aware** action of the in-basis-restricted orbit-rep +/// generator `M` at momentum `sector`, plus the `(m, s)`/`μ` selection data +/// — from ONE action pass over the basis. +/// +/// `cols[c]` holds `(row, χ_k(g_{cnt_q}) · v_q)` for every action output +/// Pauli `q` of `L*(basis[c])` whose orbit rep `r_q` is in `basis` at index +/// `row`; outputs whose rep is out of basis are dropped. This is the +/// expensive part of the orbit-rep dynamics (`compute_action_terms`, +/// [`Sector::canonicalize_phase`]). +/// +/// Unlike [`build_mf_cols`], `per_col[c].0` sums only the retained +/// in-basis entries — the exact column 1-norm of the restricted `M`, not an +/// upper bound: several distinct outputs `q` can share one rep, so the +/// out-of-basis magnitudes are not attributable to a column of `M`. `diag` +/// accumulates for the same reason. +fn build_orbit_rep_cols( + spec: &LindbladSpec, + basis: &[Word], + index: &FxHashMap, + sector: Sector<'_>, +) -> (Cols>, PerCol>) { + basis + .par_iter() + .enumerate() + .map_init( + || { + ( + Vec::::with_capacity(spec.n_qubits()), + Vec::::with_capacity(128), + FxHashMap::>::with_capacity_and_hasher( + 128, + FxBuildHasher::default(), + ), + ) + }, + |(s1, s2, lm), (c, r)| { + let terms = spec.compute_action_terms(r, s1, s2, lm); + let mut out = Vec::with_capacity(terms.len()); + let mut raw = 0.0; + let mut diag = Complex::new(0.0, 0.0); + for (q, v) in terms.iter() { + let (r_q, phase) = sector.canonicalize_phase(q); + if let Some(&row) = index.get(&r_q) { + let val = phase * *v; + raw += val.norm(); + if row as usize == c { + diag += val; + } + out.push((row, val)); + } + } + (out, (raw, diag)) + }, + ) + .unzip() +} + /// Borrowed CSC-style view of an in-basis-restricted generator `M`, backed /// by a cached per-column action computed once per expm call /// ([`build_mf_cols`]). `dot` performs the CSC matvec `y = M·x` against the cache; the @@ -214,14 +279,58 @@ where } } +/// Shared tail of every matrix-free expm: from the cached per-column action +/// derive the diagonal shift `μ = tr(M)/n` and a bound on the column 1-norm +/// of `M − μ·I` (`raw − |diag| + |diag − μ|` per column), pick the Taylor +/// partition via `select` from `‖dt·(M−μI)‖₁`, and hand everything to +/// [`quspin_expm::ExpmOp::from_parts`]. Returns `exp(dt · M) · coeffs`. +/// +/// `select` maps `‖dt·(M−μI)‖₁` to `(m*, s, backward-error tol)`; the two +/// call sites differ only in that choice. +fn expm_apply_cached( + cols: &Cols, + per_col: &PerCol, + dt: f64, + coeffs: &[T], + select: impl FnOnce(f64) -> (u32, u32, f64), +) -> Vec +where + T: ExpmComputation + + Coeff + + PartialEq + + num::Zero + + AddAssign + + Mul + + Sub + + Div + + From + + Sum, +{ + let n = cols.len(); + let trace: T = per_col.iter().map(|(_, d)| *d).sum(); + let mu = trace / n as f64; + let onenorm = per_col + .iter() + .map(|&(raw, diag)| raw - diag.mag() + (diag - mu).mag()) + .fold(0.0_f64, f64::max); + let (m_star, s, expm_tol) = select(dt.abs() * onenorm); + + let mut v = coeffs.to_vec(); + let op = CscOp { cols, dim: n }; + let expm = + quspin_expm::ExpmOp::from_parts(op, T::from(dt), mu, s as usize, m_star as usize, expm_tol); + expm.apply(ndarray::ArrayViewMut1::from(v.as_mut_slice())) + .expect("expm apply"); + v +} + /// Compute `exp(dt · M) · coeffs` for the in-basis-restricted generator /// `M`, matrix-free, via `quspin-expm`. Returns a fresh `Vec` of length /// `basis.len()`. /// -/// One matrix-free pass extracts the diagonal shift `μ = tr(M)/n` and the -/// exact column 1-norm of `M − μ·I`; from `‖dt·(M−μI)‖₁` we pick the Taylor -/// partition `(m*, s)` and hand everything to -/// [`quspin_expm::ExpmOp::from_parts`]. +/// ONE action pass builds the CSC cache `cols` (reused across every matvec) +/// and, in the same pass, the `(raw, diag)` data the `μ`/1-norm selection +/// needs; [`expm_apply_cached`] does the rest. pub(crate) fn expm_apply_mf( spec: &LindbladSpec, basis: &[Word], @@ -229,26 +338,12 @@ pub(crate) fn expm_apply_mf( coeffs: &[f64], drop_tol: f64, ) -> Vec { - let n = basis.len(); - if n == 0 { + if basis.is_empty() { return Vec::new(); } - - // ONE action pass: build the CSC cache `cols` (reused across every matvec) - // and, in the same pass, `per_col = (raw, diag)` for the `μ`/1-norm - // selection. `raw = Σ|coeff|` (all outputs), `diag` = coeff of the - // output == input term. From these: `trace = Σ diag`, `μ = trace/n`, and - // the column 1-norm of `M − μ·I` is `raw − |diag| + |diag − μ|`. let index = build_basis_index(basis); let (cols, per_col) = build_mf_cols(spec, basis, &index); - let trace: f64 = per_col.iter().map(|(_, d)| *d).sum(); - let mu = trace / n as f64; - let onenorm = per_col - .iter() - .map(|(raw, diag)| raw - diag.abs() + (diag - mu).abs()) - .fold(0.0_f64, f64::max); - // Pick the Taylor backward-error tolerance to match the basis truncation: // when the PC prunes coarsely (drop_tol >= 1e-4) a double-precision exp is // ~10 orders more accurate than the state it acts on, so the relaxed @@ -256,24 +351,41 @@ pub(crate) fn expm_apply_mf( // lower-degree Taylor polynomial and cuts the SpMV count with no effect on // the truncated result. At tight/zero drop_tol we keep double precision so // the exact-reference paths (orbit-rep / merged) still agree bit-for-bit. - let t_norm = dt.abs() * onenorm; - let (m_star, s, expm_tol) = if drop_tol >= 1e-4 { - let (m, s) = expm::select_ms_loose(t_norm); - (m, s, 1e-6_f64) - } else { - let (m, s) = expm::select_ms(t_norm); - (m, s, 1e-12_f64) - }; + expm_apply_cached(&cols, &per_col, dt, coeffs, |t_norm| { + if drop_tol >= 1e-4 { + let (m, s) = expm::select_ms_loose(t_norm); + (m, s, 1e-6) + } else { + let (m, s) = expm::select_ms(t_norm); + (m, s, 1e-12) + } + }) +} - let mut v = coeffs.to_vec(); - let op = CscOp { - cols: &cols, - dim: n, - }; - let expm = quspin_expm::ExpmOp::from_parts(op, dt, mu, s as usize, m_star as usize, expm_tol); - expm.apply(ndarray::ArrayViewMut1::from(v.as_mut_slice())) - .expect("expm apply"); - v +/// Compute `exp(dt · M) · coeffs` for the in-basis-restricted **orbit-rep** +/// generator `M` at momentum `sector`, via `quspin-expm`. Returns a fresh +/// `Vec>` of length `basis.len()`. +/// +/// The expensive phase-aware action is computed ONCE here (via +/// [`build_orbit_rep_cols`]) and reused, CSC-style, across every +/// Krylov–Taylor matvec, exactly as on the real path. +pub(crate) fn expm_apply_orbit_rep( + spec: &LindbladSpec, + basis: &[Word], + sector: Sector<'_>, + dt: f64, + coeffs: &[Complex], +) -> Vec> { + if basis.is_empty() { + return Vec::new(); + } + let index = build_basis_index(basis); + let (cols, per_col) = build_orbit_rep_cols(spec, basis, &index, sector); + + expm_apply_cached(&cols, &per_col, dt, coeffs, |t_norm| { + let (m, s) = expm::select_ms(t_norm); + (m, s, 1e-12) + }) } /// `exp(dt · M) · b` where `M` is the REAL in-basis-restricted generator but diff --git a/crates/ppvm-lindblad/src/orbit_rep.rs b/crates/ppvm-lindblad/src/orbit_rep.rs deleted file mode 100644 index c72e3da48..000000000 --- a/crates/ppvm-lindblad/src/orbit_rep.rs +++ /dev/null @@ -1,463 +0,0 @@ -// SPDX-FileCopyrightText: 2026 The PPVM Authors -// SPDX-License-Identifier: Apache-2.0 - -//! Per-step orbit-representative evolution under translation symmetry. -//! -//! The state lives entirely in **orbit-rep form** throughout: `basis` -//! contains only canonical translation-orbit representatives, and -//! `coeffs` are complex (one per rep). The dynamics `L*` is computed -//! with **phase-aware action** — for each output Pauli `q`, we -//! canonicalize `q` to its orbit rep `r_q` with shift counter `cnt_q`, -//! and accumulate `χ_k(g_{cnt_q}) · v · c_r` (where `v` is the matrix -//! element of `L*` between input rep `r` and output `q`). -//! -//! The orbit-rep basis is ~|G|× smaller than the full-basis -//! representation, throughout the entire evolution. -//! -//! The phase-aware action is genuinely **complex** (because of the -//! `χ_k(g)` phase factors). Rather than materialise a sparse matrix, the -//! per-column action — for each input rep, the list of `(row, χ_k·v)` -//! pairs for the in-basis outputs — is computed **once per expm call** -//! (via [`build_orbit_rep_cols`]) and then reused, CSC-style, across -//! every Krylov–Taylor matvec driving the external `quspin-expm` engine. -//! -//! ## Limitations -//! -//! - Caller is responsible for ensuring the input basis is in orbit-rep -//! form (i.e. each entry is the canonical representative of its -//! translation orbit). Use [`canonicalize_basis_to_rep`] if needed. -//! - The momentum sector `k_modes` is fixed for the duration of one -//! pc_step call. To compute a full site-resolved profile, call -//! `pc_step_orbit_rep` once per momentum mode and inverse-Fourier -//! the results. - -use crate::mf_expm::CscOp; -use crate::{Error, LindbladSpec}; -use fxhash::{FxBuildHasher, FxHashMap}; -use num::Complex; -use ppvm_pauli_sum::symmetry::TranslationGroup; -use quspin_expm::ExpmOp; -use rayon::prelude::*; - -// Word type re-exported from lib.rs. -use crate::Word; - -/// Replace each entry of `basis` with its canonical orbit -/// representative under `group`. Pure rewrite; coefficients are -/// untouched. Useful to enforce the orbit-rep invariant before calling -/// [`pc_step_orbit_rep`]. -/// -/// Does NOT deduplicate — if multiple input entries collapse to the -/// same rep, both are kept (caller should run a merge afterwards). -pub fn canonicalize_basis_to_rep(basis: &mut [Word], group: &TranslationGroup) { - for w in basis.iter_mut() { - *w = group.canonicalize(w); - } -} - -/// Build the per-column phase-aware action of the in-basis-restricted -/// orbit-rep generator `M` at momentum sector `k_modes`. -/// -/// Returns, for each input rep `basis[c]` (column `c`), the list of -/// `(row, χ_k(g_{cnt_q}) · v_q)` pairs for every action output Pauli `q` -/// of `L*(basis[c])` whose orbit rep `r_q` is in `basis` at index `row`. -/// Outputs not in `basis` are dropped. This is the expensive part of the -/// orbit-rep dynamics (`compute_action_terms`, `canonicalize_with_shift`, -/// `character`); it is computed once and reused by the CSC-style matvec -/// in [`CscOp`]. -pub(crate) fn build_orbit_rep_cols( - spec: &LindbladSpec, - basis: &[Word], - index: &FxHashMap, - group: &TranslationGroup, - k_modes: &[i32], -) -> Vec)>> { - basis - .par_iter() - .map_init( - || { - ( - Vec::::with_capacity(spec.n_qubits()), - Vec::::with_capacity(128), - FxHashMap::>::with_capacity_and_hasher( - 128, - FxBuildHasher::default(), - ), - ) - }, - |(s1, s2, lm), r| { - let terms = spec.compute_action_terms(r, s1, s2, lm); - let mut out = Vec::with_capacity(terms.len()); - for (q, v) in terms.iter() { - let (r_q, cnt_q) = group.canonicalize_with_shift(q); - if let Some(&row) = index.get(&r_q) { - let phase = group.character(k_modes, &cnt_q); - out.push((row, phase * *v)); - } - } - out - }, - ) - .collect() -} - -/// Compute `exp(dt · M) · coeffs` for the in-basis-restricted orbit-rep -/// generator `M` at momentum sector `k_modes`, via `quspin-expm`. Returns -/// a fresh `Vec>` of length `basis.len()`. -/// -/// The expensive phase-aware action is computed ONCE here (via -/// [`build_orbit_rep_cols`]) and reused, CSC-style, across every Krylov– -/// Taylor matvec (see [`CscOp`]). One pass over the cached columns -/// extracts the diagonal shift `μ = tr(M)/n` and a valid upper bound on -/// the column 1-norm of `M − μ·I`; from `‖dt·(M−μI)‖₁` we pick the Taylor -/// partition `(m*, s)` and hand everything to -/// [`quspin_expm::ExpmOp::from_parts`] (mirroring -/// [`crate::mf_expm::expm_apply_mf`]). -pub(crate) fn expm_apply_orbit_rep_cached( - spec: &LindbladSpec, - basis: &[Word], - group: &TranslationGroup, - k_modes: &[i32], - dt: f64, - coeffs: &[Complex], -) -> Vec> { - let n = basis.len(); - if n == 0 { - return Vec::new(); - } - - let index = crate::build_basis_index(basis); - let cols = build_orbit_rep_cols(spec, basis, &index, group, k_modes); - - // One pass for the per-column `(raw, diag)` used by the `μ`/1-norm - // selection: `raw = Σ|val|` (upper bound on the absolute column sum), - // `diag = M[c,c]`. From these: `trace = Σ diag`, `μ = trace/n`, and an - // upper bound on the column 1-norm of `M − μ·I`: `raw − |diag| + |diag − μ|`. - let per_col: Vec<(f64, Complex)> = cols - .par_iter() - .enumerate() - .map(|(c, col)| { - let mut raw = 0.0_f64; - let mut diag = Complex::new(0.0, 0.0); - for &(row, val) in col.iter() { - raw += val.norm(); - if row as usize == c { - diag += val; - } - } - (raw, diag) - }) - .collect(); - - let trace: Complex = per_col.iter().map(|(_, d)| *d).sum(); - let mu = trace / n as f64; - let onenorm = per_col - .iter() - .map(|(raw, diag)| raw - diag.norm() + (diag - mu).norm()) - .fold(0.0_f64, f64::max); - - let (m_star, s) = crate::expm::select_ms(dt.abs() * onenorm); - - let mut v = coeffs.to_vec(); - let op = CscOp { - cols: &cols, - dim: n, - }; - let e = ExpmOp::from_parts( - op, - Complex::new(dt, 0.0), - mu, - s as usize, - m_star as usize, - 1e-12_f64, - ); - e.apply(ndarray::ArrayViewMut1::from(v.as_mut_slice())) - .expect("expm apply (orbit-rep cached)"); - v -} - -/// 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 `k_modes`). -/// -/// 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, cnt_q)`. -/// 2. If `r_q` NOT in `basis` and NOT in `protected`: -/// `merged[r_q] += χ_k(g_{cnt_q}) · v_q · c_r`. -/// -/// Returns `(r_q, sum)` pairs for all candidates with nonzero sum. -/// -/// The live candidate map is capped to the *available room* -/// `room = max_basis − basis.len()` (the reps we could actually add), -/// applied during accumulation: input reps are processed in descending -/// `|c|` order and after each chunk only the `room` largest-`|sum|` -/// candidates are kept. A large `max_basis` (room ≥ all candidates) -/// disables the cap — the near-exact case. -#[allow(clippy::too_many_arguments)] -pub fn leakage_orbit_rep( - spec: &LindbladSpec, - basis: &[Word], - coeffs: &[Complex], - protected: &[Word], - group: &TranslationGroup, - k_modes: &[i32], - max_basis: usize, -) -> Result)>, Error> { - if basis.len() != coeffs.len() { - return Err(Error::LengthMismatch { - what: "basis and coeffs", - a: basis.len(), - b: coeffs.len(), - }); - } - let in_basis: FxHashMap<&Word, ()> = basis.iter().map(|w| (w, ())).collect(); - let protected_set: FxHashMap<&Word, ()> = protected.iter().map(|w| (w, ())).collect(); - - // Descending sort by |c|: process largest-magnitude contributors first - // so the running room-cap keeps the right entries. - let mut order: Vec = (0..basis.len()).collect(); - order.sort_by(|&a, &b| { - coeffs[b] - .norm() - .partial_cmp(&coeffs[a].norm()) - .unwrap_or(std::cmp::Ordering::Equal) - }); - - const CHUNK_SIZE: usize = 4096; - let room = max_basis.saturating_sub(basis.len()); - let mut merged: FxHashMap> = FxHashMap::default(); - for chunk_indices in order.chunks(CHUNK_SIZE) { - let local: Vec)>> = chunk_indices - .par_iter() - .map_init( - || { - ( - Vec::::with_capacity(spec.n_qubits()), - Vec::::with_capacity(128), - FxHashMap::>::with_capacity_and_hasher( - 128, - FxBuildHasher::default(), - ), - ) - }, - |(s1, s2, lm), &i| { - let r = &basis[i]; - let c_r = coeffs[i]; - let terms = spec.compute_action_terms(r, s1, s2, lm); - let mut out = Vec::with_capacity(terms.len()); - for (q, v) in terms.iter() { - let (r_q, cnt_q) = group.canonicalize_with_shift(q); - if !in_basis.contains_key(&r_q) && !protected_set.contains_key(&r_q) { - let phase = group.character(k_modes, &cnt_q); - out.push((r_q, phase * *v * c_r)); - } - } - out - }, - ) - .collect(); - for v in local { - for (k, val) in v { - *merged.entry(k).or_insert(Complex::new(0.0, 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 = merged.values().map(|v| v.norm()).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.norm() >= cutoff); - } - } - } - Ok(merged.into_iter().filter(|(_, c)| c.norm() > 0.0).collect()) -} - -/// Per-step orbit-rep predictor-corrector evolution. -/// -/// All state lives in orbit-rep form throughout. Each pc step does: -/// 1. Phase-aware leakage from `(basis, coeffs)`; append the largest -/// leakage reps, up to the admission room. -/// 2. Predictor: cached-action expm ([`expm_apply_orbit_rep_cached`]). -/// 3. Phase-aware leakage from the predicted state; append further reps. -/// 4. Corrector: cached-action expm from the pre-step coefficients on -/// the doubly-enlarged basis. -/// 5. Prune `|c| < drop_tol`, then trim to the top-`max_basis` reps by -/// `|c|`; protected reps never dropped. -/// -/// `max_basis` is a hard rank cap on the live orbit-rep basis: enrichment -/// adds at most `max_basis − basis.len()` of the largest leakage reps, the -/// leakage map is capped to the same room, and the post-step basis is -/// trimmed to the top-`max_basis` by `|c|`. Pass a large value (e.g. -/// `usize::MAX`) for the near-exact, uncapped case. `drop_tol` additionally -/// prunes by magnitude. -/// -/// `basis` is assumed to contain only canonical orbit representatives. -/// If not, [`canonicalize_basis_to_rep`] should be called first. -#[allow(clippy::too_many_arguments)] -pub fn pc_step_orbit_rep( - spec: &LindbladSpec, - basis: &mut Vec, - coeffs: &mut Vec>, - dt: f64, - protected: &[Word], - group: &TranslationGroup, - k_modes: &[i32], - cfg: &crate::PcStepConfig, -) -> Result<(), Error> { - let crate::PcStepConfig { - max_basis, - admit_basis, - drop_tol, - tau_add, - .. - } = *cfg; - // Admission bound, mirroring the real-space `pc_step`: enrichment may - // grow the live basis to `admit` >= `max_basis`; the final - // `cap_basis_complex` keeps the top-`max_basis` reps by evolved |coeff| - // over the whole union (rank displacement). With `admit_basis = None` - // admission is bounded by `max_basis` itself and membership turnover - // requires `drop_tol > 0`. - let admit = admit_basis.unwrap_or(max_basis).max(max_basis); - let tau_add = tau_add.unwrap_or(0.0); - // 1. First-hop phase-aware leakage. - let mut leak = leakage_orbit_rep(spec, basis, coeffs, protected, group, k_modes, admit)?; - if tau_add > 0.0 { - leak.retain(|(_, c)| c.norm() > tau_add); - } - add_leakage_capped_complex(basis, coeffs, leak, admit); - // 2. Predictor: cached-action expm (the phase-aware action is built - // once via `build_orbit_rep_cols`). - let coeffs_predict = expm_apply_orbit_rep_cached(spec, basis, group, k_modes, dt, coeffs); - // 3. Second-hop leakage from predicted state. - let mut leak2 = leakage_orbit_rep( - spec, - basis, - &coeffs_predict, - protected, - group, - k_modes, - admit, - )?; - drop(coeffs_predict); - if tau_add > 0.0 { - leak2.retain(|(_, c)| c.norm() > tau_add); - } - add_leakage_capped_complex(basis, coeffs, leak2, admit); - // 4. Corrector: cache-the-action expm from pre-step state (basis grew). - *coeffs = expm_apply_orbit_rep_cached(spec, basis, group, k_modes, dt, coeffs); - // 5. Prune by magnitude, then rank-cap to max_basis. - if drop_tol > 0.0 { - prune_basis_complex_local(basis, coeffs, drop_tol, protected); - } - cap_basis_complex(basis, coeffs, max_basis, protected); - Ok(()) -} - -/// Complex analogue of `crate::add_leakage_capped`: add the largest leakage -/// reps to the basis, up to the available room `room = max_basis − -/// basis.len()`, so the in-step orbit-rep basis never exceeds `max_basis`. -/// New reps get coefficient 0; the surrounding expm fills them. No -/// magnitude filter — the top-`room` by `|leakage|` are added. -fn add_leakage_capped_complex( - basis: &mut Vec, - coeffs: &mut Vec>, - mut leak: Vec<(Word, Complex)>, - max_basis: usize, -) { - let room = max_basis.saturating_sub(basis.len()); - if leak.len() > room { - if room > 0 { - leak.select_nth_unstable_by(room - 1, |a, b| { - b.1.norm() - .partial_cmp(&a.1.norm()) - .unwrap_or(std::cmp::Ordering::Equal) - }); - } - leak.truncate(room); - } - for (w, _) in leak { - basis.push(w); - coeffs.push(Complex::new(0.0, 0.0)); - } -} - -/// Complex analogue of `crate::cap_basis`: keep only the `max_basis` -/// largest-`|c|` reps (protected reps always kept), dropping the rest. -/// A `max_basis` large enough to cover the whole basis is a no-op. -fn cap_basis_complex( - basis: &mut Vec, - coeffs: &mut Vec>, - max_basis: usize, - protected: &[Word], -) { - if basis.len() <= max_basis { - return; - } - let protected_set: fxhash::FxHashSet<&Word> = protected.iter().collect(); - let n_prot = basis.iter().filter(|w| protected_set.contains(w)).count(); - let slots = max_basis.saturating_sub(n_prot); - let mut mags: Vec = basis - .iter() - .zip(coeffs.iter()) - .filter(|(w, _)| !protected_set.contains(w)) - .map(|(_, c)| c.norm()) - .collect(); - let cutoff = if slots == 0 { - f64::INFINITY - } else if slots >= mags.len() { - return; - } else { - let k = slots - 1; - mags.select_nth_unstable_by(k, |a, b| { - b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal) - }); - mags[k] - }; - let mut write = 0; - for read in 0..basis.len() { - if protected_set.contains(&basis[read]) || coeffs[read].norm() >= cutoff { - if write != read { - basis.swap(write, read); - coeffs.swap(write, read); - } - write += 1; - } - } - basis.truncate(write); - coeffs.truncate(write); -} - -/// Complex analogue of `crate::prune_basis`: drop reps with `|c| < -/// drop_tol`, never dropping `protected` reps. No-op when `drop_tol <= 0`. -fn prune_basis_complex_local( - basis: &mut Vec, - coeffs: &mut Vec>, - drop_tol: f64, - protected: &[Word], -) { - if drop_tol <= 0.0 { - return; - } - let protected_set: fxhash::FxHashSet<&Word> = protected.iter().collect(); - let mut write = 0; - for read in 0..basis.len() { - if coeffs[read].norm() >= drop_tol || protected_set.contains(&basis[read]) { - if write != read { - basis.swap(write, read); - coeffs.swap(write, read); - } - write += 1; - } - } - basis.truncate(write); - coeffs.truncate(write); -} diff --git a/crates/ppvm-lindblad/src/scalar.rs b/crates/ppvm-lindblad/src/scalar.rs new file mode 100644 index 000000000..9ae3f4175 --- /dev/null +++ b/crates/ppvm-lindblad/src/scalar.rs @@ -0,0 +1,43 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! The coefficient scalar of a Pauli-sum basis: real (`f64`) on the +//! plain adaptive path, complex on the momentum-sector orbit-rep path. +//! +//! Every truncation and 1-norm decision in this crate only ever needs a +//! magnitude and a zero, so the two paths share one implementation +//! parameterised by [`Coeff`] instead of a real and a complex copy. + +use num::Complex; + +/// A Pauli-sum coefficient: `f64` or `Complex`. +pub(crate) trait Coeff: Copy + Send + Sync { + /// Absolute value (`f64::abs`) / modulus (`Complex::norm`). + fn mag(self) -> f64; + + fn zero() -> Self; +} + +impl Coeff for f64 { + #[inline] + fn mag(self) -> f64 { + self.abs() + } + + #[inline] + fn zero() -> Self { + 0.0 + } +} + +impl Coeff for Complex { + #[inline] + fn mag(self) -> f64 { + self.norm() + } + + #[inline] + fn zero() -> Self { + Complex::new(0.0, 0.0) + } +} diff --git a/crates/ppvm-lindblad/src/sector.rs b/crates/ppvm-lindblad/src/sector.rs new file mode 100644 index 000000000..ab588f7aa --- /dev/null +++ b/crates/ppvm-lindblad/src/sector.rs @@ -0,0 +1,75 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Momentum sectors of a translation group, and the phase-aware +//! canonicalization that drives orbit-representative evolution. +//! +//! On the orbit-rep path the state lives entirely in **orbit-rep form** +//! throughout: the basis contains only canonical translation-orbit +//! representatives and the coefficients are complex (one per rep). The +//! dynamics `L*` is computed with **phase-aware action** — for each +//! output Pauli `q`, we canonicalize `q` to its orbit rep `r_q` with +//! shift counter `cnt_q`, and accumulate `χ_k(g_{cnt_q}) · v · c_r` +//! (where `v` is the matrix element of `L*` between input rep `r` and +//! output `q`). [`Sector::canonicalize_phase`] is that step. +//! +//! The orbit-rep basis is ~`|G|`× smaller than the full-basis +//! representation, throughout the entire evolution. +//! +//! ## Limitations +//! +//! - Callers are responsible for ensuring the input basis is in +//! orbit-rep form (i.e. each entry is the canonical representative of +//! its translation orbit). Use [`canonicalize_basis_to_rep`] if +//! needed. +//! - A [`Sector`] is fixed for the duration of one +//! [`LindbladSpec::pc_step_orbit_rep`](crate::LindbladSpec::pc_step_orbit_rep) +//! call. To compute a full site-resolved profile, call it once per +//! momentum mode and inverse-Fourier the results. + +use crate::Word; +use num::Complex; +use ppvm_pauli_sum::symmetry::TranslationGroup; + +/// A momentum sector of a translation group: the group `G` together with +/// one integer mode index per generator. The wavenumber along generator +/// `g` is `2π · k_modes[g] / group.generator_order(g)`; `k_modes = [0, …]` +/// is the trivial sector. +/// +/// The two halves are meaningless apart — every phase-aware routine +/// needs both — so they travel as one value. +#[derive(Clone, Copy)] +pub struct Sector<'a> { + group: &'a TranslationGroup, + k_modes: &'a [i32], +} + +impl<'a> Sector<'a> { + pub fn new(group: &'a TranslationGroup, k_modes: &'a [i32]) -> Self { + Self { group, k_modes } + } + + /// Canonicalize `q` to its orbit representative `r_q` and return it + /// alongside the character phase `χ_k(g_{cnt_q})` of the group + /// element that maps `q` to `r_q`. The phase weights the matrix + /// element of `L*` when it is accumulated onto `r_q`. + #[inline] + pub fn canonicalize_phase(&self, q: &Word) -> (Word, Complex) { + let (rep, counter) = self.group.canonicalize_with_shift(q); + let phase = self.group.character(self.k_modes, &counter); + (rep, phase) + } +} + +/// Replace each entry of `basis` with its canonical orbit +/// representative under `group`. Pure rewrite; coefficients are +/// untouched. Useful to enforce the orbit-rep invariant before calling +/// [`LindbladSpec::pc_step_orbit_rep`](crate::LindbladSpec::pc_step_orbit_rep). +/// +/// Does NOT deduplicate — if multiple input entries collapse to the +/// same rep, both are kept (caller should run a merge afterwards). +pub fn canonicalize_basis_to_rep(basis: &mut [Word], group: &TranslationGroup) { + for w in basis.iter_mut() { + *w = group.canonicalize(w); + } +} diff --git a/crates/ppvm-lindblad/src/step.rs b/crates/ppvm-lindblad/src/step.rs index fe77ce83f..d0e380a0e 100644 --- a/crates/ppvm-lindblad/src/step.rs +++ b/crates/ppvm-lindblad/src/step.rs @@ -3,10 +3,12 @@ //! Predictor-corrector adaptive step `O ← exp(dt·L*) O`. +use crate::sector::Sector; use crate::spec::LindbladSpec; +use crate::truncate::{add_leakage_capped, cap_basis, prune_basis}; use crate::word::Word; use crate::{Error, PcStepConfig, mf_expm}; -use fxhash::FxHashSet; +use num::Complex; use std::time::Instant; /// Per-phase timing breakdown (microseconds) returned by @@ -48,99 +50,6 @@ impl Phase { } } -/// Compact `basis` / `coeffs` in place: drop entries whose absolute -/// coefficient is below `drop_tol` unless the word appears in `protected`. -/// No-op when `drop_tol ≤ 0`. -fn prune_basis(basis: &mut Vec, coeffs: &mut Vec, drop_tol: f64, protected: &[Word]) { - if drop_tol <= 0.0 { - return; - } - debug_assert_eq!(basis.len(), coeffs.len()); - let protected_set: FxHashSet<&Word> = protected.iter().collect(); - let mut write = 0; - for read in 0..basis.len() { - if coeffs[read].abs() >= drop_tol || protected_set.contains(&basis[read]) { - if write != read { - basis.swap(write, read); - coeffs.swap(write, read); - } - write += 1; - } - } - basis.truncate(write); - coeffs.truncate(write); -} - -/// Global max-basis cap (PauliStrings.jl-style top-M trim): keep only the -/// `max_basis` largest-|coeff| terms (protected strings always kept), -/// dropping the rest. Rank-based total-basis bound; dual of `drop_tol`. -/// A `max_basis` large enough to cover the whole basis is a no-op. -fn cap_basis(basis: &mut Vec, coeffs: &mut Vec, max_basis: usize, protected: &[Word]) { - if basis.len() <= max_basis { - return; - } - let protected_set: FxHashSet<&Word> = protected.iter().collect(); - let n_prot = basis.iter().filter(|w| protected_set.contains(w)).count(); - let slots = max_basis.saturating_sub(n_prot); - let mut mags: Vec = basis - .iter() - .zip(coeffs.iter()) - .filter(|(w, _)| !protected_set.contains(w)) - .map(|(_, c)| c.abs()) - .collect(); - let cutoff = if slots == 0 { - f64::INFINITY - } else if slots >= mags.len() { - return; - } else { - let k = slots - 1; - mags.select_nth_unstable_by(k, |a, b| { - b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal) - }); - mags[k] - }; - let mut write = 0; - for read in 0..basis.len() { - if protected_set.contains(&basis[read]) || coeffs[read].abs() >= cutoff { - if write != read { - basis.swap(write, read); - coeffs.swap(write, read); - } - write += 1; - } - } - basis.truncate(write); - coeffs.truncate(write); -} - -/// Add the largest leakage strings to the basis, up to the available room -/// `room = max_basis − basis.len()` — so the in-step basis (hence the -/// expm/leakage peak memory) never exceeds `max_basis`. New strings get -/// coefficient 0; the surrounding expm fills them. No magnitude filter: the -/// top-`room` by `|leakage|` are added (a large `max_basis` adds them all). -fn add_leakage_capped( - basis: &mut Vec, - coeffs: &mut Vec, - mut leak: Vec<(Word, f64)>, - max_basis: usize, -) { - let room = max_basis.saturating_sub(basis.len()); - if leak.len() > room { - if room > 0 { - leak.select_nth_unstable_by(room - 1, |a, b| { - b.1.abs() - .partial_cmp(&a.1.abs()) - .unwrap_or(std::cmp::Ordering::Equal) - }); - } - leak.truncate(room); - } - for (w, _) in leak { - basis.push(w); - coeffs.push(0.0); - } -} - impl LindbladSpec { /// One predictor-corrector step `O ← exp(dt·L*) O` in the adaptive /// real-coefficient Pauli basis: first-hop leakage admission, predictor @@ -269,4 +178,76 @@ impl LindbladSpec { fn expm_step(&self, basis: &[Word], dt: f64, b: &[f64], drop_tol: f64) -> Vec { mf_expm::expm_apply_mf(self, basis, dt, b, drop_tol) } + + /// One predictor-corrector step in **orbit-rep form** at momentum + /// `sector`: the same five phases as [`Self::pc_step`], but the basis + /// holds only canonical translation-orbit representatives, the + /// coefficients are complex, and the `L*` action is phase-aware (see + /// [`crate::sector`]). The basis stays ~`|G|`× smaller than the + /// equivalent full-basis complex evolution, every step. + /// + /// `max_basis` is a hard rank cap on the live orbit-rep basis: + /// enrichment adds at most `admit − basis.len()` of the largest + /// leakage reps, the leakage map is capped to the same room, and the + /// post-step basis is trimmed to the top-`max_basis` reps by `|c|`. + /// Pass a large value (e.g. `usize::MAX`) for the near-exact, + /// uncapped case. `drop_tol` additionally prunes by magnitude. + /// `protected` reps are never dropped. + /// + /// `basis` is assumed to contain only canonical orbit + /// representatives. If not, call + /// [`canonicalize_basis_to_rep`](crate::canonicalize_basis_to_rep) + /// first. + pub fn pc_step_orbit_rep( + &self, + basis: &mut Vec, + coeffs: &mut Vec>, + dt: f64, + protected: &[Word], + sector: Sector<'_>, + cfg: &PcStepConfig, + ) -> Result<(), Error> { + let PcStepConfig { + max_basis, + admit_basis, + drop_tol, + tau_add, + .. + } = *cfg; + // Admission bound, mirroring `pc_step_inner`: enrichment may grow + // the live basis to `admit` >= `max_basis`; the final `cap_basis` + // keeps the top-`max_basis` reps by evolved |coeff| over the whole + // union (rank displacement). With `admit_basis = None` admission is + // bounded by `max_basis` itself and membership turnover requires + // `drop_tol > 0`. + let admit = admit_basis.unwrap_or(max_basis).max(max_basis); + let tau_add = tau_add.unwrap_or(0.0); + + // 1. First-hop phase-aware leakage. + let mut leak = self.leakage_orbit_rep(basis, coeffs, protected, sector, admit)?; + if tau_add > 0.0 { + leak.retain(|(_, c)| c.norm() > tau_add); + } + add_leakage_capped(basis, coeffs, leak, admit); + + // 2. Predictor: the phase-aware action is built once and reused + // across every matvec. + let coeffs_predict = mf_expm::expm_apply_orbit_rep(self, basis, sector, dt, coeffs); + + // 3. Second-hop leakage from the predicted state. + let mut leak2 = self.leakage_orbit_rep(basis, &coeffs_predict, protected, sector, admit)?; + drop(coeffs_predict); + if tau_add > 0.0 { + leak2.retain(|(_, c)| c.norm() > tau_add); + } + add_leakage_capped(basis, coeffs, leak2, admit); + + // 4. Corrector: redo from the pre-step state (the basis grew). + *coeffs = mf_expm::expm_apply_orbit_rep(self, basis, sector, dt, coeffs); + + // 5. Prune by magnitude, then rank-cap to max_basis. + prune_basis(basis, coeffs, drop_tol, protected); + cap_basis(basis, coeffs, max_basis, protected); + Ok(()) + } } diff --git a/crates/ppvm-lindblad/src/tests.rs b/crates/ppvm-lindblad/src/tests.rs index 3ad4899d5..5687c0b31 100644 --- a/crates/ppvm-lindblad/src/tests.rs +++ b/crates/ppvm-lindblad/src/tests.rs @@ -155,15 +155,14 @@ fn pc_step_orbit_rep_matches_full_basis_projection() { let mut cr = coeffs_full.clone(); canonicalize_pauli_sum_complex(&mut br, &mut cr, &group, &k); // Evolve in orbit-rep form (max_basis large ⇒ full enrichment). + let sector = Sector::new(&group, &k); for _ in 0..n_steps { - orbit_rep::pc_step_orbit_rep( - &spec, + spec.pc_step_orbit_rep( &mut br, &mut cr, dt, &protected, - &group, - &k, + sector, &PcStepConfig { max_basis: 10_000_000, ..Default::default() diff --git a/crates/ppvm-lindblad/src/truncate.rs b/crates/ppvm-lindblad/src/truncate.rs new file mode 100644 index 000000000..c96a927fb --- /dev/null +++ b/crates/ppvm-lindblad/src/truncate.rs @@ -0,0 +1,153 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Basis truncation and enrichment: the magnitude prune, the rank cap, +//! and the capped leakage admission shared by every `pc_step` variant. +//! +//! All three are generic over the coefficient scalar ([`Coeff`]), so the +//! real adaptive path and the complex orbit-rep path run the same code. + +use crate::scalar::Coeff; +use crate::word::Word; +use fxhash::{FxHashMap, FxHashSet}; + +/// Basis indices in descending coefficient magnitude. Leakage +/// accumulation walks the basis in this order so the running room-cap +/// keeps the entries most likely to be the true largest contributors. +pub(crate) fn order_by_desc_mag(coeffs: &[T]) -> Vec { + let mut order: Vec = (0..coeffs.len()).collect(); + order.sort_by(|&a, &b| desc_by_mag(coeffs[a], coeffs[b])); + order +} + +/// Keep only the `room` largest-magnitude entries of a live leakage +/// candidate map — `room` being the number of strings we could actually +/// admit to the basis, so there is no point tracking more. Applied after +/// each accumulation chunk. +pub(crate) fn cap_map_to_room(merged: &mut FxHashMap, room: usize) { + if merged.len() <= room { + return; + } + if room == 0 { + merged.clear(); + return; + } + let mut mags: Vec = merged.values().map(|v| v.mag()).collect(); + let k = room.min(mags.len() - 1); + let cutoff = nth_largest(&mut mags, k); + merged.retain(|_, v| v.mag() >= cutoff); +} + +/// Compact `basis` / `coeffs` in place: drop entries whose coefficient +/// magnitude is below `drop_tol` unless the word appears in `protected`. +/// No-op when `drop_tol ≤ 0`. +pub(crate) fn prune_basis( + basis: &mut Vec, + coeffs: &mut Vec, + drop_tol: f64, + protected: &[Word], +) { + if drop_tol <= 0.0 { + return; + } + debug_assert_eq!(basis.len(), coeffs.len()); + let protected_set: FxHashSet<&Word> = protected.iter().collect(); + retain_in_place(basis, coeffs, |w, c| { + c.mag() >= drop_tol || protected_set.contains(w) + }); +} + +/// Global max-basis cap (PauliStrings.jl-style top-M trim): keep only the +/// `max_basis` largest-magnitude terms (protected strings always kept), +/// dropping the rest. Rank-based total-basis bound; dual of `drop_tol`. +/// A `max_basis` large enough to cover the whole basis is a no-op. +pub(crate) fn cap_basis( + basis: &mut Vec, + coeffs: &mut Vec, + max_basis: usize, + protected: &[Word], +) { + if basis.len() <= max_basis { + return; + } + let protected_set: FxHashSet<&Word> = protected.iter().collect(); + let n_prot = basis.iter().filter(|w| protected_set.contains(w)).count(); + let slots = max_basis.saturating_sub(n_prot); + let mut mags: Vec = basis + .iter() + .zip(coeffs.iter()) + .filter(|(w, _)| !protected_set.contains(w)) + .map(|(_, c)| c.mag()) + .collect(); + let cutoff = if slots == 0 { + f64::INFINITY + } else if slots >= mags.len() { + return; + } else { + nth_largest(&mut mags, slots - 1) + }; + retain_in_place(basis, coeffs, |w, c| { + protected_set.contains(w) || c.mag() >= cutoff + }); +} + +/// Add the largest leakage strings to the basis, up to the available room +/// `room = max_basis − basis.len()` — so the in-step basis (hence the +/// expm/leakage peak memory) never exceeds `max_basis`. New strings get +/// coefficient 0; the surrounding expm fills them. No magnitude filter: the +/// top-`room` by `|leakage|` are added (a large `max_basis` adds them all). +pub(crate) fn add_leakage_capped( + basis: &mut Vec, + coeffs: &mut Vec, + mut leak: Vec<(Word, T)>, + max_basis: usize, +) { + let room = max_basis.saturating_sub(basis.len()); + if leak.len() > room { + if room > 0 { + leak.select_nth_unstable_by(room - 1, |a, b| desc_by_mag(a.1, b.1)); + } + leak.truncate(room); + } + for (w, _) in leak { + basis.push(w); + coeffs.push(T::zero()); + } +} + +/// Keep the `basis`/`coeffs` entries satisfying `keep`, preserving order, +/// by swapping survivors down and truncating. +fn retain_in_place( + basis: &mut Vec, + coeffs: &mut Vec, + mut keep: impl FnMut(&Word, &T) -> bool, +) { + let mut write = 0; + for read in 0..basis.len() { + if keep(&basis[read], &coeffs[read]) { + if write != read { + basis.swap(write, read); + coeffs.swap(write, read); + } + write += 1; + } + } + basis.truncate(write); + coeffs.truncate(write); +} + +/// The `k`-th largest element of `mags` (0-indexed), via a partial sort. +/// Reorders `mags`. Panics if `k >= mags.len()`. +fn nth_largest(mags: &mut [f64], k: usize) -> f64 { + mags.select_nth_unstable_by(k, |a, b| { + b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal) + }); + mags[k] +} + +/// Descending comparison by magnitude, NaN-tolerant. +fn desc_by_mag(a: T, b: T) -> std::cmp::Ordering { + b.mag() + .partial_cmp(&a.mag()) + .unwrap_or(std::cmp::Ordering::Equal) +} diff --git a/crates/ppvm-python-native/src/lindblad.rs b/crates/ppvm-python-native/src/lindblad.rs index 7c7206201..88799a321 100644 --- a/crates/ppvm-python-native/src/lindblad.rs +++ b/crates/ppvm-python-native/src/lindblad.rs @@ -391,7 +391,7 @@ impl LindbladSpec { tau_add: Option, ) -> PyResult> { use num::Complex; - use ppvm_lindblad::orbit_rep; + use ppvm_lindblad::{Sector, canonicalize_basis_to_rep}; let n_q = self.inner.n_qubits(); let basis_view = basis.as_array(); @@ -422,25 +422,24 @@ impl LindbladSpec { ))); } if canonicalize_first { - orbit_rep::canonicalize_basis_to_rep(&mut basis_words, group.core()); + canonicalize_basis_to_rep(&mut basis_words, group.core()); } - orbit_rep::pc_step_orbit_rep( - &self.inner, - &mut basis_words, - &mut coeffs_vec, - dt, - &protected_words, - group.core(), - k_slice, - &ppvm_lindblad::PcStepConfig { - max_basis, - admit_basis, - drop_tol, - tau_add, - num_threads: None, - }, - ) - .map_err(map_err)?; + self.inner + .pc_step_orbit_rep( + &mut basis_words, + &mut coeffs_vec, + dt, + &protected_words, + Sector::new(group.core(), k_slice), + &ppvm_lindblad::PcStepConfig { + max_basis, + admit_basis, + drop_tol, + tau_add, + num_threads: None, + }, + ) + .map_err(map_err)?; let m = basis_words.len(); let mut out_basis = vec![0u8; m * n_q]; From 4fa2dbd09afaff1861504af5614019ca29ab26da Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 2 Sep 2026 11:52:38 +0200 Subject: [PATCH 08/13] refactor(pauli-sum): lift momentum_merge out of the PyO3 macro MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `PauliSum.momentum_merge` carried ~90 lines of numerics inside a `macro_rules!` body: gather two real sums into a `HashMap`, fold onto orbit reps, rescale, split back into real/imaginary parts. That made it untestable from Rust (only reachable through Python) and expanded it once per `Config` variant. The algorithm now lives in `ppvm_pauli_sum::symmetry` as `momentum_merge_pauli_sum_pair`, beside `canonicalize_pauli_sum_complex` which it reuses, and parallel to `symmetry_merge_pauli_sum`. The PyO3 method keeps only its qubit-count / momentum-length validation (so the Python error messages are unchanged) and one call. New Rust tests pin the two things the wrapper used to hide: - `momentum_merge_pair_matches_symmetry_merge_at_k0`: on free orbits the k=0 pair merge agrees with the independent real-coefficient `symmetry_merge_pauli_sum` path — this is what pins the `|G|` rescale. - `momentum_merge_pair_matches_summing_projector`: for every k on a 4-chain, the rep coefficient equals `Σ_orbit χ_k(g) · c`, computed from the group API. - two `should_panic` tests for the qubit-count and momentum-length asserts. No behavior change. `cargo clippy --workspace --all-targets` is clean and the Python suite still passes. Co-Authored-By: Claude Opus 5 (1M context) --- crates/ppvm-pauli-sum/src/symmetry/mod.rs | 5 +- .../ppvm-pauli-sum/src/symmetry/momentum.rs | 89 ++++++++++- crates/ppvm-pauli-sum/src/symmetry/tests.rs | 147 ++++++++++++++++++ crates/ppvm-python-native/src/interface.rs | 54 +------ 4 files changed, 245 insertions(+), 50 deletions(-) diff --git a/crates/ppvm-pauli-sum/src/symmetry/mod.rs b/crates/ppvm-pauli-sum/src/symmetry/mod.rs index 30d6d4c01..99ada0665 100644 --- a/crates/ppvm-pauli-sum/src/symmetry/mod.rs +++ b/crates/ppvm-pauli-sum/src/symmetry/mod.rs @@ -68,7 +68,10 @@ mod momentum; pub use group::TranslationGroup; pub use merge::{canonicalize_pauli_sum, symmetry_merge_pauli_sum}; -pub use momentum::{SectorCheckError, canonicalize_pauli_sum_complex, check_momentum_sector}; +pub use momentum::{ + SectorCheckError, canonicalize_pauli_sum_complex, check_momentum_sector, + momentum_merge_pauli_sum_pair, +}; #[cfg(test)] mod tests; diff --git a/crates/ppvm-pauli-sum/src/symmetry/momentum.rs b/crates/ppvm-pauli-sum/src/symmetry/momentum.rs index ad29facdc..952d459a2 100644 --- a/crates/ppvm-pauli-sum/src/symmetry/momentum.rs +++ b/crates/ppvm-pauli-sum/src/symmetry/momentum.rs @@ -1,10 +1,11 @@ // SPDX-FileCopyrightText: 2026 The PPVM Authors // SPDX-License-Identifier: Apache-2.0 +use crate::sum::PauliSum; use fxhash::{FxHashMap, FxHashSet}; use num::Complex; use ppvm_pauli_word::word::PauliWord; -use ppvm_traits::{HashFinalize, PauliStorage}; +use ppvm_traits::{ACMapAddAssign, ACMapBase, ACMapIter, Config, HashFinalize, PauliStorage}; use std::f64::consts::PI; use std::hash::BuildHasher; @@ -143,6 +144,92 @@ pub fn canonicalize_pauli_sum_complex( } } +/// Momentum-sector merge of a complex operator carried as a **real +/// pair**: `re` and `im` are the real and imaginary parts of +/// `O = re + i·im`. Both are overwritten in place with the +/// orbit-representative form of `O` projected onto momentum sector +/// `k_modes`. +/// +/// This is the momentum-sector counterpart of +/// [`super::symmetry_merge_pauli_sum`], and generalizes it to `k ≠ 0` +/// while keeping real coefficients on both sums — the only complex +/// arithmetic is the internal character-weighted fold, which reuses +/// [`canonicalize_pauli_sum_complex`]. +/// +/// [`canonicalize_pauli_sum_complex`] carries a `1/|orbit|` prefactor +/// (it *averages* over orbit members); we rescale by `group.order()` so +/// that the result is the *summing* projector, matching +/// [`super::symmetry_merge_pauli_sum`]. +/// +/// Entries whose rescaled component is exactly zero are dropped, so a +/// purely real operator leaves `im` empty. +/// +/// # Panics +/// +/// If `re` and `im` disagree on qubit count, if either disagrees with +/// `group.n_qubits()`, or if `k_modes.len() != group.n_generators()`. +pub fn momentum_merge_pauli_sum_pair( + re: &mut PauliSum, + im: &mut PauliSum, + group: &TranslationGroup, + k_modes: &[i32], +) where + T: Config, Coeff = f64>, + T::Map: ACMapAddAssign>, + for<'a> T::Map: ACMapIter<'a, Item = (&'a PauliWord, &'a f64)>, + A: PauliStorage, + S: BuildHasher + Clone + Default + HashFinalize, +{ + assert_eq!( + re.n_qubits(), + im.n_qubits(), + "real and imaginary parts disagree on qubit count" + ); + assert_eq!( + re.n_qubits(), + group.n_qubits(), + "PauliSum qubit count {} != group qubit count {}", + re.n_qubits(), + group.n_qubits() + ); + assert_eq!( + k_modes.len(), + group.n_generators(), + "k_modes length {} != number of generators {}", + k_modes.len(), + group.n_generators() + ); + + // Gather both real components into `word -> re + i·im`. + let mut combined: FxHashMap, Complex> = FxHashMap::default(); + for (word, coeff) in re.data().iter() { + combined.entry(*word).or_insert(Complex::new(0.0, 0.0)).re += *coeff; + } + for (word, coeff) in im.data().iter() { + combined.entry(*word).or_insert(Complex::new(0.0, 0.0)).im += *coeff; + } + let mut basis = Vec::with_capacity(combined.len()); + let mut coeffs = Vec::with_capacity(combined.len()); + for (word, coeff) in combined { + basis.push(word); + coeffs.push(coeff); + } + + canonicalize_pauli_sum_complex(&mut basis, &mut coeffs, group, k_modes); + + let scale = group.order() as f64; + re.data_mut().clear(); + im.data_mut().clear(); + for (word, coeff) in basis.into_iter().zip(coeffs) { + if coeff.re != 0.0 { + *re += (word, coeff.re * scale); + } + if coeff.im != 0.0 { + *im += (word, coeff.im * scale); + } + } +} + /// Verify that a `(basis, complex_coeffs)` Pauli sum lies entirely in /// the momentum sector `k_modes` under `group`. /// diff --git a/crates/ppvm-pauli-sum/src/symmetry/tests.rs b/crates/ppvm-pauli-sum/src/symmetry/tests.rs index 64c1deec0..ee9981e60 100644 --- a/crates/ppvm-pauli-sum/src/symmetry/tests.rs +++ b/crates/ppvm-pauli-sum/src/symmetry/tests.rs @@ -2,6 +2,7 @@ // SPDX-License-Identifier: Apache-2.0 use super::*; +use crate::sum::PauliSum; use fxhash::FxHashMap; use num::Complex; use ppvm_pauli_word::word::PauliWord; @@ -452,6 +453,152 @@ fn pauli_sum_symmetry_merge_matches_plain_trotter() { ); } +/// Build the `(re, im)` real pair of the momentum-`k` eigenstate +/// `O_k = Σ_a e^{-2πi k a / n} Z_a` on an `n`-site chain. +fn seed_z_momentum_pair(n: usize, k: i32) -> (PauliSum, PauliSum) +where + Cfg: ppvm_traits::Config, + PauliSum: for<'s> std::ops::AddAssign<(&'s str, f64)>, +{ + let mut re: PauliSum = PauliSum::builder().n_qubits(n).build(); + let mut im: PauliSum = PauliSum::builder().n_qubits(n).build(); + for a in 0..n { + let mut s: Vec = vec!['I'; n]; + s[a] = 'Z'; + let st: String = s.into_iter().collect(); + let phase = -2.0 * PI * (k as f64) * (a as f64) / (n as f64); + re += (st.as_str(), phase.cos()); + im += (st.as_str(), phase.sin()); + } + (re, im) +} + +/// At `k = 0` on free orbits, the phase-aware pair merge must agree +/// with the independent real-coefficient `symmetry_merge_pauli_sum` +/// code path — this is what pins the `|G|` rescale (the momentum +/// projector averages; `symmetry_merge` sums). +#[test] +fn momentum_merge_pair_matches_symmetry_merge_at_k0() { + use crate::config::indexmap::ByteFxHashF64; + use crate::prelude::*; + + type Cfg = ByteFxHashF64<1>; + let n = 4usize; + let group = TranslationGroup::chain_1d(n); + + // Σ_j Z_j and Σ_j X_j X_{j+1}: both free orbits (|orbit| = |G|). + let mut reference: PauliSum = PauliSum::builder().n_qubits(n).build(); + let mut re: PauliSum = PauliSum::builder().n_qubits(n).build(); + let mut im: PauliSum = PauliSum::builder().n_qubits(n).build(); + for j in 0..n { + let mut z: Vec = vec!['I'; n]; + z[j] = 'Z'; + let mut xx: Vec = vec!['I'; n]; + xx[j] = 'X'; + xx[(j + 1) % n] = 'X'; + for (s, c) in [(z, 1.5), (xx, -0.25)] { + let st: String = s.into_iter().collect(); + reference += (st.as_str(), c); + re += (st.as_str(), c); + } + } + // `im` needs an entry to exist; a zero coefficient must not survive. + im += ("IIII", 0.0); + + symmetry_merge_pauli_sum(&mut reference, &group); + momentum_merge_pauli_sum_pair(&mut re, &mut im, &group, &[0]); + + let expected: FxHashMap<_, f64> = reference.iter().map(|(w, c)| (*w, *c)).collect(); + let got: FxHashMap<_, f64> = re.iter().map(|(w, c)| (*w, *c)).collect(); + assert_eq!(expected.len(), got.len(), "basis sizes differ"); + for (w, &c) in &expected { + let g = *got.get(w).unwrap_or_else(|| panic!("missing rep {w:?}")); + assert!( + (c - g).abs() < 1e-12, + "rep {w:?}: symmetry_merge gave {c}, momentum_merge gave {g}" + ); + } + assert_eq!(im.len(), 0, "a purely real input must leave `im` empty"); +} + +/// A momentum-`k` eigenstate folds onto a single orbit rep whose +/// coefficient is the *summing* projector `Σ_{p ∈ orbit} χ_k(g_p) · c_p`, +/// computed here directly from the group API. +#[test] +fn momentum_merge_pair_matches_summing_projector() { + use crate::config::indexmap::ByteFxHashF64; + + type Cfg = ByteFxHashF64<1>; + let n = 4usize; + let group = TranslationGroup::chain_1d(n); + + for k in 0..n as i32 { + let (mut re, mut im) = seed_z_momentum_pair::(n, k); + // Coefficient of each orbit member before merging. + let before: FxHashMap> = (0..n) + .map(|a| { + let mut s: Vec = vec!['I'; n]; + s[a] = 'Z'; + let phase = -2.0 * PI * (k as f64) * (a as f64) / (n as f64); + ( + word(&s.into_iter().collect::()), + Complex::from_polar(1.0, phase), + ) + }) + .collect(); + + momentum_merge_pauli_sum_pair(&mut re, &mut im, &group, &[k]); + + let got_re: FxHashMap<_, f64> = re.iter().map(|(w, c)| (*w, *c)).collect(); + let got_im: FxHashMap<_, f64> = im.iter().map(|(w, c)| (*w, *c)).collect(); + assert!( + !got_re.is_empty() || !got_im.is_empty(), + "k={k}: merged away" + ); + + let rep = group.canonicalize(&word(&{ + let mut s: Vec = vec!['I'; n]; + s[0] = 'Z'; + s.into_iter().collect::() + })); + // Σ over the orbit of χ_k(g) · c_{g·rep}. + let mut expected = Complex::new(0.0, 0.0); + for (member, counter) in group.orbit_with_counters(&rep) { + expected += group.character(&[k], &counter) * before[&member]; + } + let got = Complex::new( + got_re.get(&rep).copied().unwrap_or(0.0), + got_im.get(&rep).copied().unwrap_or(0.0), + ); + assert!( + (got - expected).norm() < 1e-12, + "k={k}: rep {rep:?} expected {expected:?}, got {got:?}" + ); + } +} + +#[test] +#[should_panic(expected = "k_modes length 2 != number of generators 1")] +fn momentum_merge_pair_rejects_wrong_momentum_length() { + use crate::config::indexmap::ByteFxHashF64; + + type Cfg = ByteFxHashF64<1>; + let group = TranslationGroup::chain_1d(4); + let (mut re, mut im) = seed_z_momentum_pair::(4, 0); + momentum_merge_pauli_sum_pair(&mut re, &mut im, &group, &[0, 0]); +} + +#[test] +#[should_panic(expected = "PauliSum qubit count 4 != group qubit count 3")] +fn momentum_merge_pair_rejects_qubit_count_mismatch() { + use crate::config::indexmap::ByteFxHashF64; + + type Cfg = ByteFxHashF64<1>; + let group = TranslationGroup::chain_1d(3); + let (mut re, mut im) = seed_z_momentum_pair::(4, 0); + momentum_merge_pauli_sum_pair(&mut re, &mut im, &group, &[0]); +} + #[test] #[should_panic(expected = "generator 0 order must be nonzero")] fn rejects_zero_generator_order() { diff --git a/crates/ppvm-python-native/src/interface.rs b/crates/ppvm-python-native/src/interface.rs index 54813f130..e84c29fe5 100644 --- a/crates/ppvm-python-native/src/interface.rs +++ b/crates/ppvm-python-native/src/interface.rs @@ -105,15 +105,15 @@ macro_rules! create_interface_symmetry_methods { group: &crate::symmetry::TranslationGroup, momentum: Vec, ) -> pyo3::PyResult<()> { - let n_g = group.core().n_qubits(); + let n_q = group.core().n_qubits(); for (label, n) in [ ("self", self.inner.n_qubits()), ("other", other.inner.n_qubits()), ] { - if n != n_g { + if n != n_q { return Err(pyo3::exceptions::PyValueError::new_err(format!( "{label} PauliSum has {n} qubits but the \ - TranslationGroup acts on {n_g}", + TranslationGroup acts on {n_q}", ))); } } @@ -124,54 +124,12 @@ macro_rules! create_interface_symmetry_methods { group.core().n_generators(), ))); } - // Gather both real components into word -> (re + i·im). - let mut combined: std::collections::HashMap< - <$type as Config>::PauliWordType, - num::Complex, - > = std::collections::HashMap::new(); - for (w, v) in self.inner.data().iter() { - combined - .entry(w.clone()) - .or_insert(num::Complex::new(0.0, 0.0)) - .re += *v; - } - for (w, v) in other.inner.data().iter() { - combined - .entry(w.clone()) - .or_insert(num::Complex::new(0.0, 0.0)) - .im += *v; - } - let mut basis = Vec::with_capacity(combined.len()); - let mut coeffs = Vec::with_capacity(combined.len()); - for (w, c) in combined { - basis.push(w); - coeffs.push(c); - } - // Character-weighted fold onto orbit reps. - // `canonicalize_pauli_sum_complex` carries a 1/|G| prefactor; - // we rescale by |G| so the merge is the *summing* projector - // (like `symmetry_merge`): idempotent on already-merged input, - // hence stable under merging after every Trotter step. - ppvm_pauli_sum::symmetry::canonicalize_pauli_sum_complex( - &mut basis, - &mut coeffs, + ppvm_pauli_sum::symmetry::momentum_merge_pauli_sum_pair( + &mut self.inner, + &mut other.inner, group.core(), &momentum, ); - let scale = group.core().order() as f64; - // Write the real/imag parts back into the two sums. - self.inner.data_mut().clear(); - other.inner.data_mut().clear(); - for (w, c) in basis.into_iter().zip(coeffs.into_iter()) { - let re = c.re * scale; - let im = c.im * scale; - if re != 0.0 { - self.inner += (w.clone(), re); - } - if im != 0.0 { - other.inner += (w, im); - } - } Ok(()) } } From 4841ffc9ac5d1943e53e0ca1305fe0148234225d Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 2 Sep 2026 11:53:48 +0200 Subject: [PATCH 09/13] test: cover the untested Python symmetry surface MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Five public Python entry points shipped with no Python test at all; `momentum_merge` was the only one covered. Adds: - `test_symmetry_merge.py` — `TranslationGroup` constructors, properties, `from_generators` validation, `canonicalize`, and `PauliSum.symmetry_merge` (orbit summing, distinct orbits, idempotency, coefficient-sum conservation, qubit-count mismatch). - `test_symmetry_arrays.py` — the three `_core` array functions (`canonicalize_basis_arr`, `canonicalize_basis_arr_complex`, `check_momentum_sector_arr`) against numpy references derived from the group action, plus every validation path. `check_momentum_sector_arr` was previously exported and referenced from docstrings but never exercised. - `test/lindblad/test_pc_step_orbit_rep.py` — the `pc_step_orbit_rep` binding: orbit-rep evolution vs a dense full-space numpy matrix exponential projected at the end (k = 0, 1, 2), `canonicalize_first`, the `max_basis` rank cap, `protected` reps, input validation, and returned dtypes/shapes. Also adds an xfail (strict) documenting a real bug found while writing these: `momentum_merge` rescales by `|G|`, but the projector averages over the `|orbit|` DISTINCT members, so an orbit with a non-trivial stabilizer is amplified by `|G|/|orbit|` on every merge — `ZZZZ` on a 4-chain grows 4x per merge, `ZIZI` 2x. The existing idempotency test only seeds free orbits, where the two factors coincide, which is why it passes. Co-Authored-By: Claude Opus 5 (1M context) --- .../test/lindblad/test_pc_step_orbit_rep.py | 222 ++++++++++++++++++ ppvm-python/test/test_momentum_merge.py | 26 ++ ppvm-python/test/test_symmetry_arrays.py | 218 +++++++++++++++++ ppvm-python/test/test_symmetry_merge.py | 142 +++++++++++ 4 files changed, 608 insertions(+) create mode 100644 ppvm-python/test/lindblad/test_pc_step_orbit_rep.py create mode 100644 ppvm-python/test/test_symmetry_arrays.py create mode 100644 ppvm-python/test/test_symmetry_merge.py diff --git a/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py b/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py new file mode 100644 index 000000000..fdb234cdb --- /dev/null +++ b/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py @@ -0,0 +1,222 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 + +"""Orbit-representative predictor-corrector evolution through the Python +binding (:meth:`Lindbladian.pc_step_orbit_rep`). + +The state lives entirely in orbit-rep form: the basis holds only canonical +translation-orbit representatives and the coefficients are complex. For a +translation-invariant Lindbladian, evolving in orbit-rep form and projecting +onto the momentum sector at the *end* of a full-space evolution must agree — +that is the content of the projection theorem, and it is what the first test +checks against a dense numpy matrix exponential, independent of the Rust +Al-Mohy & Higham implementation. +""" + +from __future__ import annotations + +import cmath + +import numpy as np +import pytest + +from ppvm import Lindbladian +from ppvm._core import TranslationGroup, canonicalize_basis_arr_complex + +from ._helpers import all_strings + +_CODE = {"I": 0, "X": 1, "Z": 2, "Y": 3} +_CHAR = {v: k for k, v in _CODE.items()} + + +def basis_arr(strings, n): + arr = np.zeros((len(strings), n), dtype=np.uint8) + for i, s in enumerate(strings): + arr[i] = [_CODE[c] for c in s] + return arr + + +def string(row): + return "".join(_CHAR[int(c)] for c in row) + + +def to_dict(basis, coeffs): + return {string(w): c for w, c in zip(basis, coeffs, strict=True)} + + +def momentum(*modes): + return np.array(modes, dtype=np.int32) + + +def xy_chain_pbc(n, gamma): + """Translation-invariant XY chain with PBC plus uniform Z dephasing.""" + h_terms = [] + for j in range(n): + nxt = (j + 1) % n + for op in "XY": + s = ["I"] * n + s[j] = op + s[nxt] = op + h_terms.append(("".join(s), 1.0)) + jumps = [("I" * j + "Z" + "I" * (n - j - 1), gamma) for j in range(n)] + return Lindbladian(n, h_terms, jumps) + + +def z_momentum_seed(n, k): + """``O_k = Σ_a e^{-2πi k a / n} Z_a`` as ``(basis_arr, complex coeffs)``.""" + words = ["I" * a + "Z" + "I" * (n - a - 1) for a in range(n)] + coeffs = np.array([cmath.exp(-2j * cmath.pi * k * a / n) for a in range(n)]) + return basis_arr(words, n), coeffs + + +def _dense_expm(A, terms=40): + """``exp(A)`` for a real matrix, by Taylor series with scaling and squaring. + + The full-space Lindbladian is not diagonalizable in general (numpy's + ``eig`` returns a singular eigenvector matrix here), and the test suite + deliberately has no scipy dependency — so the reference is this direct + series, independent of the Rust Al-Mohy & Higham implementation. + """ + norm = np.abs(A).sum(axis=0).max() + s = max(0, int(np.ceil(np.log2(norm))) + 1) if norm > 0 else 0 + B = A / 2**s + total = np.eye(A.shape[0]) + term = np.eye(A.shape[0]) + for k in range(1, terms + 1): + term = term @ B / k + total = total + term + for _ in range(s): + total = total @ total + return total + + +@pytest.mark.parametrize("k", [0, 1, 2]) +def test_orbit_rep_matches_dense_full_space_then_project(k): + """Orbit-rep evolution == full-space evolution projected at the end. + + Full space is all 4^n Pauli strings (n=3 -> 64), exponentiated densely + with numpy. The orbit-rep side runs with a huge ``max_basis`` so its rank + cap never binds and the only remaining difference would be a bug in the + phase-aware action. + """ + n = 3 + dt = 0.02 + n_steps = 3 + op = xy_chain_pbc(n, gamma=0.3) + group = TranslationGroup.chain_1d(n) + k_arr = momentum(k) + + # --- dense full-space reference --- + full = all_strings(n) + generator = np.zeros((len(full), len(full)), dtype=float) + rows, cols, vals = op.generator(full) + generator[rows, cols] = vals + seed_basis, seed_coeffs = z_momentum_seed(n, k) + index = {s: i for i, s in enumerate(full)} + v = np.zeros(len(full), dtype=complex) + for w, c in zip(seed_basis, seed_coeffs, strict=True): + v[index[string(w)]] = c + # `generator` is real, so exp(dt·G) is too: evolve the real and imaginary + # parts of the coefficient vector separately. + step = _dense_expm(dt * generator) + v_re, v_im = v.real.copy(), v.imag.copy() + for _ in range(n_steps): + v_re = step @ v_re + v_im = step @ v_im + v = v_re + 1j * v_im + expected = to_dict(*canonicalize_basis_arr_complex(basis_arr(full, n), v, group, k_arr)) + + # --- orbit-rep evolution --- + rep_basis, rep_coeffs = canonicalize_basis_arr_complex(seed_basis, seed_coeffs, group, k_arr) + for _ in range(n_steps): + rep_basis, rep_coeffs = op.pc_step_orbit_rep( + rep_basis, rep_coeffs, dt, 10_000_000, group, k_arr, drop_tol=0.0 + ) + got = to_dict(rep_basis, rep_coeffs) + + # Compare on the union; the dense side keeps exact zeros the orbit-rep + # side never admits, so only nonzero entries must match. + for word in set(expected) | set(got): + e = expected.get(word, 0.0) + g = got.get(word, 0.0) + assert abs(e - g) < 1e-9, f"k={k}: rep {word} dense {e} vs orbit-rep {g}" + assert any(abs(c) > 1e-6 for c in got.values()), "orbit-rep state decayed away" + + +def test_canonicalize_first_accepts_non_canonical_input(): + """The same physical state seeded on a non-canonical orbit member gives + the same evolution once ``canonicalize_first=True`` normalizes it.""" + n = 3 + dt = 0.02 + op = xy_chain_pbc(n, gamma=0.3) + group = TranslationGroup.chain_1d(n) + k_arr = momentum(0) + + seed_basis, seed_coeffs = z_momentum_seed(n, 0) + canonical, coeffs = canonicalize_basis_arr_complex(seed_basis, seed_coeffs, group, k_arr) + + ref_basis, ref_coeffs = op.pc_step_orbit_rep(canonical, coeffs, dt, 10_000_000, group, k_arr) + # Feed a shifted (non-canonical) representative of the same orbit. + shifted = np.array([[_CODE[c] for c in "IZI"]], dtype=np.uint8) + got_basis, got_coeffs = op.pc_step_orbit_rep( + shifted, coeffs, dt, 10_000_000, group, k_arr, canonicalize_first=True + ) + ref = to_dict(ref_basis, ref_coeffs) + got = to_dict(got_basis, got_coeffs) + assert ref.keys() == got.keys() + for w in ref: + assert abs(ref[w] - got[w]) < 1e-12 + + +def test_max_basis_caps_the_live_basis(): + n = 4 + op = xy_chain_pbc(n, gamma=0.1) + group = TranslationGroup.chain_1d(n) + k_arr = momentum(0) + seed_basis, seed_coeffs = z_momentum_seed(n, 0) + basis, coeffs = canonicalize_basis_arr_complex(seed_basis, seed_coeffs, group, k_arr) + for _ in range(4): + basis, coeffs = op.pc_step_orbit_rep(basis, coeffs, 0.05, 6, group, k_arr) + assert basis.shape[0] <= 6 + assert basis.shape == (len(coeffs), n) + + +def test_protected_reps_are_never_dropped(): + n = 4 + op = xy_chain_pbc(n, gamma=0.1) + group = TranslationGroup.chain_1d(n) + k_arr = momentum(0) + seed_basis, seed_coeffs = z_momentum_seed(n, 0) + basis, coeffs = canonicalize_basis_arr_complex(seed_basis, seed_coeffs, group, k_arr) + protected = basis.copy() + keep = {string(w) for w in protected} + for _ in range(3): + # max_basis=1 with a huge drop_tol would wipe everything unprotected. + basis, coeffs = op.pc_step_orbit_rep( + basis, coeffs, 0.05, 1, group, k_arr, drop_tol=1e3, protected_arr=protected + ) + assert keep <= {string(w) for w in basis} + + +def test_pc_step_orbit_rep_validates_inputs(): + n = 3 + op = xy_chain_pbc(n, gamma=0.0) + group = TranslationGroup.chain_1d(n) + basis, coeffs = z_momentum_seed(n, 0) + with pytest.raises(ValueError, match="momentum has 2 entries but group has 1 generators"): + op.pc_step_orbit_rep(basis, coeffs, 0.01, 100, group, momentum(0, 0)) + with pytest.raises(ValueError, match="coeffs has length 2 but basis has 3 rows"): + op.pc_step_orbit_rep(basis, coeffs[:2], 0.01, 100, group, momentum(0)) + + +def test_returns_complex_arrays_of_matching_shape(): + n = 3 + op = xy_chain_pbc(n, gamma=0.2) + group = TranslationGroup.chain_1d(n) + k_arr = momentum(1) + seed_basis, seed_coeffs = z_momentum_seed(n, 1) + basis, coeffs = canonicalize_basis_arr_complex(seed_basis, seed_coeffs, group, k_arr) + out_basis, out_coeffs = op.pc_step_orbit_rep(basis, coeffs, 0.01, 500, group, k_arr) + assert out_basis.dtype == np.uint8 + assert out_coeffs.dtype == np.complex128 + assert out_basis.shape == (len(out_coeffs), n) diff --git a/ppvm-python/test/test_momentum_merge.py b/ppvm-python/test/test_momentum_merge.py index 0546deb11..6dccfc84d 100644 --- a/ppvm-python/test/test_momentum_merge.py +++ b/ppvm-python/test/test_momentum_merge.py @@ -91,6 +91,32 @@ def test_momentum_merge_idempotent(k): assert max(abs(once.get(x, 0j) - twice.get(x, 0j)) for x in keys) < 1e-12 +@pytest.mark.xfail( + reason="momentum_merge rescales by |G| but the projector averages over the " + "|orbit| DISTINCT members, so orbits with a non-trivial stabilizer are " + "amplified by |G|/|orbit| on every merge. `_seed_pair` only produces free " + "orbits, which is why test_momentum_merge_idempotent passes.", + strict=True, +) +@pytest.mark.parametrize("word, orbit_size", [("ZZZZ", 1), ("ZIZI", 2)]) +def test_momentum_merge_idempotent_on_stabilized_orbit(word, orbit_size): + """Idempotency must hold for every orbit, not just the free ones. + + ``ZZZZ`` is translation-invariant (orbit size 1) and ``ZIZI`` has period 2, + so on a 4-site chain they pick up factors of 4 and 2 per merge. + """ + n = 4 + g = TranslationGroup.chain_1d(n) + PA = PauliSum.new(n, [(word, 1.0)], min_abs_coeff=0.0, max_pauli_weight=n) + PB = PauliSum.new(n, [(word, 0.0)], min_abs_coeff=0.0, max_pauli_weight=n) + PA.momentum_merge(PB, g, [0]) + once = _to_complex_dict(PA, PB) + PA.momentum_merge(PB, g, [0]) + twice = _to_complex_dict(PA, PB) + keys = set(once) | set(twice) + assert max(abs(once.get(x, 0j) - twice.get(x, 0j)) for x in keys) < 1e-12 + + def test_momentum_merge_projects_out_other_sectors(): """Merging a pure sector-k operator in sector k' != k gives ~zero.""" n = 4 diff --git a/ppvm-python/test/test_symmetry_arrays.py b/ppvm-python/test/test_symmetry_arrays.py new file mode 100644 index 000000000..11ec2aa59 --- /dev/null +++ b/ppvm-python/test/test_symmetry_arrays.py @@ -0,0 +1,218 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 + +"""Tests for the array-form symmetry primitives on the ``(basis_arr, coeffs)`` +representation used by ``Lindbladian.pc_step_arr``: + +- ``canonicalize_basis_arr`` — plain real merge (sums colliding coefficients) +- ``canonicalize_basis_arr_complex`` — momentum-sector projection (averages + over the distinct orbit members with the character weight) +- ``check_momentum_sector_arr`` — validation that an input really lies in the + sector it is about to be projected onto + +References are computed here in numpy from the group action, independent of +the Rust merge routines. +""" + +import cmath + +import numpy as np +import pytest + +from ppvm._core import ( + TranslationGroup, + canonicalize_basis_arr, + canonicalize_basis_arr_complex, + check_momentum_sector_arr, +) + +_CODE = {"I": 0, "X": 1, "Z": 2, "Y": 3} +_CHAR = {v: k for k, v in _CODE.items()} + + +def basis_arr(strings): + return np.array([[_CODE[c] for c in s] for s in strings], dtype=np.uint8) + + +def string(row): + return "".join(_CHAR[int(c)] for c in row) + + +def to_dict(pair): + words, coeffs = pair + return {string(w): c for w, c in zip(words, coeffs, strict=True)} + + +def rep_of(group, s): + return string(group.canonicalize(np.array([_CODE[c] for c in s], dtype=np.uint8))) + + +def momentum(*modes): + """The ``_core`` free functions take momentum as an int32 array; numpy's + default integer dtype is int64, which they reject.""" + return np.array(modes, dtype=np.int32) + + +def z_strings(n): + return ["I" * j + "Z" + "I" * (n - j - 1) for j in range(n)] + + +# ── canonicalize_basis_arr (real, k=0) ─────────────────────────────────────── +def test_canonicalize_basis_arr_sums_collisions(): + n = 4 + g = TranslationGroup.chain_1d(n) + words = z_strings(n) + coeffs = np.array([1.0, 2.0, 3.0, 4.0]) + merged = to_dict(canonicalize_basis_arr(basis_arr(words), coeffs, g)) + assert merged == pytest.approx({rep_of(g, words[0]): 10.0}) + + +def test_canonicalize_basis_arr_matches_manual_grouping(): + """Reference: group rows by their rep in numpy and sum.""" + n = 4 + g = TranslationGroup.chain_1d(n) + rng = np.random.default_rng(3) + words = [*z_strings(n), "XXII", "IXXI", "IIXX", "XIIX", "ZZZZ"] + coeffs = rng.normal(size=len(words)) + + expected: dict[str, float] = {} + for w, c in zip(words, coeffs, strict=True): + expected[rep_of(g, w)] = expected.get(rep_of(g, w), 0.0) + c + + merged = to_dict(canonicalize_basis_arr(basis_arr(words), coeffs, g)) + assert merged.keys() == expected.keys() + for w in expected: + assert merged[w] == pytest.approx(expected[w]) + + +def test_canonicalize_basis_arr_validates_shapes(): + g = TranslationGroup.chain_1d(4) + with pytest.raises(ValueError, match="3 qubits per row but group acts on 4"): + canonicalize_basis_arr(basis_arr(["ZII"]), np.array([1.0]), g) + with pytest.raises(ValueError, match="coeffs has length 2 but basis has 1 rows"): + canonicalize_basis_arr(basis_arr(["ZIII"]), np.array([1.0, 2.0]), g) + + +# ── canonicalize_basis_arr_complex (momentum sectors) ──────────────────────── +def _momentum_seed(n, k): + """``O_k = Σ_a e^{-2πi k a / n} Z_a`` as ``(basis_arr, coeffs)``.""" + words = z_strings(n) + coeffs = np.array([cmath.exp(-2j * cmath.pi * k * a / n) for a in range(n)]) + return basis_arr(words), coeffs + + +@pytest.mark.parametrize("k", [0, 1, 2, 3]) +def test_complex_merge_of_momentum_eigenstate_has_unit_rep_coefficient(k): + """The projection *averages* over the orbit, so a normalized momentum + eigenstate folds to a rep coefficient of modulus 1.""" + n = 4 + g = TranslationGroup.chain_1d(n) + words, coeffs = _momentum_seed(n, k) + merged = to_dict(canonicalize_basis_arr_complex(words, coeffs, g, momentum(k))) + assert len(merged) == 1 + assert abs(next(iter(merged.values()))) == pytest.approx(1.0) + + +def test_complex_merge_projects_out_other_sectors(): + """A pure k=1 state has zero component in every other sector.""" + n = 4 + g = TranslationGroup.chain_1d(n) + words, coeffs = _momentum_seed(n, 1) + for k_other in [0, 2, 3]: + merged = to_dict(canonicalize_basis_arr_complex(words, coeffs, g, momentum(k_other))) + for c in merged.values(): + assert abs(c) < 1e-12, f"k=1 state leaked into sector {k_other}: {c}" + + +def test_complex_merge_matches_character_average(): + """Reference: (1/|orbit|) Σ_{p in orbit} χ_k(g_p) · c_p, computed here + by walking the cyclic shifts explicitly.""" + n = 4 + k = 1 + g = TranslationGroup.chain_1d(n) + rng = np.random.default_rng(11) + words = z_strings(n) + coeffs = rng.normal(size=n) + 1j * rng.normal(size=n) + + # Z_a is the shift of Z_0 by `a`, so the character weight is e^{2πika/n}. + by_word = dict(zip(words, coeffs, strict=True)) + rep = rep_of(g, words[0]) + shift_of_rep = words.index(rep) + expected = ( + sum( + cmath.exp(2j * cmath.pi * k * ((a - shift_of_rep) % n) / n) * by_word[words[a]] + for a in range(n) + ) + / n + ) + + merged = to_dict(canonicalize_basis_arr_complex(basis_arr(words), coeffs, g, momentum(k))) + assert merged[rep] == pytest.approx(expected) + + +def test_complex_merge_validates_shapes(): + g = TranslationGroup.chain_1d(4) + words, coeffs = _momentum_seed(4, 1) + with pytest.raises(ValueError, match="momentum has 2 entries but group has 1 generators"): + canonicalize_basis_arr_complex(words, coeffs, g, momentum(0, 0)) + with pytest.raises(ValueError, match="coeffs has length 2 but basis has 4 rows"): + canonicalize_basis_arr_complex(words, coeffs[:2], g, momentum(1)) + with pytest.raises(ValueError, match="3 qubits per row but group acts on 4"): + canonicalize_basis_arr_complex(basis_arr(["ZII"]), np.array([1 + 0j]), g, momentum(0)) + + +# ── check_momentum_sector_arr ──────────────────────────────────────────────── +@pytest.mark.parametrize("k", [0, 1, 2, 3]) +def test_check_momentum_sector_accepts_eigenstate(k): + n = 4 + g = TranslationGroup.chain_1d(n) + words, coeffs = _momentum_seed(n, k) + assert check_momentum_sector_arr(words, coeffs, g, momentum(k)) is None + + +def test_check_momentum_sector_rejects_wrong_sector(): + n = 4 + g = TranslationGroup.chain_1d(n) + words, coeffs = _momentum_seed(n, 1) + with pytest.raises(ValueError, match="not in target momentum sector"): + check_momentum_sector_arr(words, coeffs, g, momentum(0)) + + +def test_check_momentum_sector_rejects_incomplete_orbit(): + """Orbit members missing from the basis count as zero, so a lone Z_0 is + not a momentum eigenstate.""" + g = TranslationGroup.chain_1d(4) + with pytest.raises(ValueError, match="not in target momentum sector"): + check_momentum_sector_arr(basis_arr(["ZIII"]), np.array([1 + 0j]), g, momentum(0)) + + +def test_check_momentum_sector_flags_incompatible_stabilizer(): + """``ZIZI`` has a period-2 stabilizer, which cannot carry k=1.""" + g = TranslationGroup.chain_1d(4) + with pytest.raises(ValueError, match="stabilizer incompatible with momentum sector"): + check_momentum_sector_arr( + basis_arr(["ZIZI", "IZIZ"]), + np.array([1 + 0j, -1 + 0j]), + g, + momentum(1), + ) + + +def test_check_momentum_sector_tolerance_is_configurable(): + n = 4 + g = TranslationGroup.chain_1d(n) + words, coeffs = _momentum_seed(n, 1) + perturbed = coeffs.copy() + perturbed[0] += 1e-6 + with pytest.raises(ValueError, match="not in target momentum sector"): + check_momentum_sector_arr(words, perturbed, g, momentum(1), 1e-9) + # Same input passes once the tolerance exceeds the perturbation. + assert check_momentum_sector_arr(words, perturbed, g, momentum(1), 1e-4) is None + + +def test_check_momentum_sector_rejects_invalid_tolerance(): + n = 4 + g = TranslationGroup.chain_1d(n) + words, coeffs = _momentum_seed(n, 0) + with pytest.raises(ValueError, match="invalid tolerance"): + check_momentum_sector_arr(words, coeffs, g, momentum(0), -1.0) diff --git a/ppvm-python/test/test_symmetry_merge.py b/ppvm-python/test/test_symmetry_merge.py new file mode 100644 index 000000000..d151cc168 --- /dev/null +++ b/ppvm-python/test/test_symmetry_merge.py @@ -0,0 +1,142 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 + +"""Tests for the ``TranslationGroup`` binding and ``PauliSum.symmetry_merge``. + +``symmetry_merge`` is the plain real-coefficient (``k=0``) merge: every Pauli +word is replaced by its canonical translation-orbit representative and +coefficients of colliding words are summed. See ``test_momentum_merge.py`` +for the phase-aware (``k != 0``) counterpart. +""" + +import numpy as np +import pytest + +from ppvm import PauliSum +from ppvm._core import TranslationGroup + +_CODE = {"I": 0, "X": 1, "Z": 2, "Y": 3} +_CHAR = {v: k for k, v in _CODE.items()} + + +def codes(s): + return np.array([_CODE[c] for c in s], dtype=np.uint8) + + +def string(arr): + return "".join(_CHAR[int(c)] for c in arr) + + +def psum(n, terms): + return PauliSum.new(n, terms, min_abs_coeff=0.0, max_pauli_weight=n) + + +# ── TranslationGroup constructors and properties ───────────────────────────── +@pytest.mark.parametrize( + "group, n_qubits, n_generators, order", + [ + (TranslationGroup.chain_1d(6), 6, 1, 6), + (TranslationGroup.torus_2d(3, 2), 6, 2, 6), + (TranslationGroup.torus_3d(2, 2, 2), 8, 3, 8), + (TranslationGroup.ladder(3, 2), 6, 1, 3), + ], +) +def test_group_shapes(group, n_qubits, n_generators, order): + assert group.n_qubits == n_qubits + assert group.n_generators == n_generators + assert group.order == order + + +def test_from_generators_matches_chain_1d(): + n = 4 + shift = [(i + 1) % n for i in range(n)] + g = TranslationGroup.from_generators(n, [shift], [n]) + ref = TranslationGroup.chain_1d(n) + assert (g.n_qubits, g.n_generators, g.order) == (ref.n_qubits, ref.n_generators, ref.order) + for s in ["ZIII", "IZII", "XYII", "IXYI"]: + assert string(g.canonicalize(codes(s))) == string(ref.canonicalize(codes(s))) + + +@pytest.mark.parametrize( + "perms, orders, message", + [ + ([[1, 0, 2, 3]], [4, 4], "same length"), + ([[1, 0, 2]], [2], "permutation length"), + ([[1, 0, 2, 9]], [2], "out of range"), + ([[1, 1, 2, 3]], [2], "duplicate target"), + ], +) +def test_from_generators_validates(perms, orders, message): + with pytest.raises(ValueError, match=message): + TranslationGroup.from_generators(4, perms, orders) + + +def test_canonicalize_is_orbit_invariant(): + g = TranslationGroup.chain_1d(4) + shifts = ["IIXY", "IXYI", "XYII", "YIIX"] + reps = {string(g.canonicalize(codes(s))) for s in shifts} + assert len(reps) == 1, "all cyclic shifts must share one representative" + # The rep is itself a member of the orbit (lex-min is over the internal + # (xbits, zbits) ordering, which isn't observable from Python). + assert reps.pop() in shifts + + +def test_canonicalize_rejects_wrong_length(): + g = TranslationGroup.chain_1d(4) + with pytest.raises(ValueError, match="length 3 but group expects 4"): + g.canonicalize(codes("IXY")) + + +# ── PauliSum.symmetry_merge ────────────────────────────────────────────────── +def test_symmetry_merge_sums_one_orbit(): + """Σ_j Z_j on a 4-chain is a single free orbit: 4 entries -> 1 with c=4.""" + n = 4 + g = TranslationGroup.chain_1d(n) + p = psum(n, [("I" * j + "Z" + "I" * (n - j - 1), 1.0) for j in range(n)]) + assert len(p.terms) == n + p.symmetry_merge(g) + assert len(p.terms) == 1 + (word, coeff) = p.terms[0] + assert coeff == pytest.approx(4.0) + assert string(g.canonicalize(codes(word))) == word + + +def test_symmetry_merge_keeps_distinct_orbits_and_weights(): + n = 4 + g = TranslationGroup.chain_1d(n) + terms = [("I" * j + "Z" + "I" * (n - j - 1), 1.0) for j in range(n)] + terms += [("I" * j + "X" + "I" * (n - j - 1), 0.25) for j in range(n)] + p = psum(n, terms) + p.symmetry_merge(g) + coeffs = sorted(c for _, c in p.terms) + assert coeffs == pytest.approx([1.0, 4.0]) + + +def test_symmetry_merge_is_idempotent(): + """A merged sum is already in orbit-rep form, so re-merging is a no-op.""" + n = 4 + g = TranslationGroup.chain_1d(n) + p = psum(n, [("I" * j + "Z" + "I" * (n - j - 1), 1.0) for j in range(n)]) + p.symmetry_merge(g) + once = sorted(p.terms) + p.symmetry_merge(g) + assert sorted(p.terms) == once + + +def test_symmetry_merge_preserves_translation_invariant_trace(): + """Merging conserves Σ_p c_p, hence any orbit-summed observable.""" + n = 4 + g = TranslationGroup.chain_1d(n) + rng = np.random.default_rng(7) + words = ["ZIII", "IZII", "IIZI", "IIIZ", "XXII", "IXXI", "IIXX", "XIIX"] + coeffs = rng.normal(size=len(words)) + p = psum(n, list(zip(words, coeffs, strict=True))) + total = sum(c for _, c in p.terms) + p.symmetry_merge(g) + assert sum(c for _, c in p.terms) == pytest.approx(total) + + +def test_symmetry_merge_rejects_qubit_count_mismatch(): + p = psum(4, [("ZIII", 1.0)]) + with pytest.raises(ValueError, match="4 qubits but the TranslationGroup acts on 3"): + p.symmetry_merge(TranslationGroup.chain_1d(3)) From 54f79beb7fb870467d0369e4fb37e2988600cf91 Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 2 Sep 2026 13:52:13 +0200 Subject: [PATCH 10/13] fix(pauli-sum): make momentum_merge the summing projector on every orbit MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `momentum_merge` rescaled the orbit-*averaged* projection by a global `group.order()`, but `canonicalize_pauli_sum_complex` divides by the number of DISTINCT orbit members. The two agree only for free orbits, so any orbit with a non-trivial stabilizer was amplified by `|G|/|orbit|` on every merge. Since the documented workflow is "merge after every Trotter step", the error compounded geometrically — on a 4-site chain `ZZZZ` grew 4x per merge and `ZIZI` 2x, so `4^steps` for a translation-invariant word: ZIII (|orbit|=4): 1.0 -> 1.0 -> 1.0 (was already correct) ZIZI (|orbit|=2): 1.0 -> 2.0 -> 4.0 now 1.0 -> 1.0 -> 1.0 ZZZZ (|orbit|=1): 1.0 -> 4.0 -> 16.0 now 1.0 -> 1.0 -> 1.0 The character-weighted fold moves into a shared `project_onto_reps`, which returns the un-normalized sum together with `|orbit|`. The two conventions are now explicit at the two call sites: `canonicalize_pauli_sum_complex` divides to average (unchanged behavior and signature), `momentum_merge_pauli_sum_pair` takes the sum as-is. Summing is what makes the merge idempotent for every orbit and makes it reduce exactly to `symmetry_merge_pauli_sum` at k=0 — which the strengthened `momentum_merge_pair_matches_symmetry_merge_at_k0` now asserts over stabilized orbits too, not just free ones. The strict xfail added in 4841ffc9 becomes a plain passing test, and every pre-existing test still passes untouched — including the exact-diagonalization Trotter checks in test_momentum_merge.py, which only ever exercised free orbits. Also corrects the docs that described the old `1/|G|` relationship. Co-Authored-By: Claude Opus 5 (1M context) --- crates/ppvm-pauli-sum/src/symmetry/mod.rs | 20 +++-- .../ppvm-pauli-sum/src/symmetry/momentum.rs | 90 ++++++++++++------- crates/ppvm-pauli-sum/src/symmetry/tests.rs | 17 ++-- crates/ppvm-python-native/src/interface.rs | 5 +- crates/ppvm-python-native/src/symmetry.rs | 5 +- ppvm-python/src/ppvm/paulisum.py | 7 +- ppvm-python/test/test_momentum_merge.py | 21 +++-- 7 files changed, 104 insertions(+), 61 deletions(-) diff --git a/crates/ppvm-pauli-sum/src/symmetry/mod.rs b/crates/ppvm-pauli-sum/src/symmetry/mod.rs index 99ada0665..22fdbb291 100644 --- a/crates/ppvm-pauli-sum/src/symmetry/mod.rs +++ b/crates/ppvm-pauli-sum/src/symmetry/mod.rs @@ -17,13 +17,19 @@ //! (`k=0`) symmetry sector, e.g. sums of single-Z operators over the //! lattice. //! -//! **Non-trivial momentum sectors (`k ≠ 0`)** are handled by -//! [`canonicalize_pauli_sum_complex`], which folds with the character -//! phase `χ_k(g)` of each translation. On the Python side, an operator in -//! sector `k` is carried as a *real pair* (real + imaginary components, two -//! real `PauliSum`s) and merged via `PauliSum.momentum_merge`, which reuses -//! this routine — letting gate-based Trotter evolution stay symmetry- -//! compressed in any momentum sector with real coefficients throughout. +//! **Non-trivial momentum sectors (`k ≠ 0`)** fold with the character +//! phase `χ_k(g)` of each translation. Two conventions share one core: +//! [`canonicalize_pauli_sum_complex`] *averages* over each orbit's +//! distinct members (`1/|orbit|`), while +//! [`momentum_merge_pauli_sum_pair`] *sums* — the convention that is +//! idempotent on every orbit and reduces exactly to +//! [`symmetry_merge_pauli_sum`] at `k = 0`. Note `|orbit| = |G|` only for +//! free orbits, so the two are **not** related by a global `|G|` factor. +//! On the Python side, an operator in sector `k` is carried as a *real +//! pair* (real + imaginary components, two real `PauliSum`s) and merged +//! via `PauliSum.momentum_merge` — letting gate-based Trotter evolution +//! stay symmetry-compressed in any momentum sector with real +//! coefficients throughout. //! //! ## Data model //! diff --git a/crates/ppvm-pauli-sum/src/symmetry/momentum.rs b/crates/ppvm-pauli-sum/src/symmetry/momentum.rs index 952d459a2..19563f27e 100644 --- a/crates/ppvm-pauli-sum/src/symmetry/momentum.rs +++ b/crates/ppvm-pauli-sum/src/symmetry/momentum.rs @@ -100,6 +100,43 @@ pub fn canonicalize_pauli_sum_complex( for (word, &coeff) in basis.iter().zip(coeffs.iter()) { *input.entry(*word).or_insert(Complex::new(0.0, 0.0)) += coeff; } + let projected = project_onto_reps(&input, group, k_modes); + basis.clear(); + coeffs.clear(); + basis.reserve(projected.len()); + coeffs.reserve(projected.len()); + for (word, (sum, orbit_size)) in projected { + basis.push(word); + coeffs.push(sum / orbit_size as f64); + } +} + +/// Character-weighted fold of `input` onto translation-orbit +/// representatives, the shared core of the two momentum-projection +/// conventions. +/// +/// Returns `rep → (Σ_{p ∈ orbit} χ_k(g_p) · c_p, |orbit|)`: the +/// **summing** projector, paired with the number of *distinct* orbit +/// members. Callers pick their convention — +/// [`canonicalize_pauli_sum_complex`] divides by `|orbit|` to average, +/// [`momentum_merge_pauli_sum_pair`] takes the sum as-is. +/// +/// `|orbit|` is `group.order()` only for free orbits; an orbit with a +/// non-trivial stabilizer has fewer distinct members, which is exactly +/// why the two conventions must not be related by a global `|G|` factor. +/// +/// Orbits whose stabilizer is incompatible with `k_modes` (the same +/// orbit member reached with different character numerators) project to +/// zero and are omitted from the output. +fn project_onto_reps( + input: &FxHashMap, Complex>, + group: &TranslationGroup, + k_modes: &[i32], +) -> FxHashMap, (Complex, usize)> +where + A: PauliStorage, + S: BuildHasher + Clone + Default + HashFinalize, +{ let reps: FxHashSet<_> = input.keys().map(|word| group.canonicalize(word)).collect(); let mut projected = FxHashMap::default(); @@ -123,25 +160,18 @@ pub fn canonicalize_pauli_sum_complex( if !compatible { continue; } - let orbit_size = members.len() as f64; - let mut rep_coeff = Complex::new(0.0, 0.0); + let orbit_size = members.len(); + let mut sum = Complex::new(0.0, 0.0); for (member, (counter, _)) in members { let coeff = input .get(&member) .copied() .unwrap_or(Complex::new(0.0, 0.0)); - rep_coeff += group.character(k_modes, &counter) * coeff / orbit_size; + sum += group.character(k_modes, &counter) * coeff; } - projected.insert(rep, rep_coeff); - } - basis.clear(); - coeffs.clear(); - basis.reserve(projected.len()); - coeffs.reserve(projected.len()); - for (w, c) in projected { - basis.push(w); - coeffs.push(c); + projected.insert(rep, (sum, orbit_size)); } + projected } /// Momentum-sector merge of a complex operator carried as a **real @@ -156,13 +186,16 @@ pub fn canonicalize_pauli_sum_complex( /// arithmetic is the internal character-weighted fold, which reuses /// [`canonicalize_pauli_sum_complex`]. /// -/// [`canonicalize_pauli_sum_complex`] carries a `1/|orbit|` prefactor -/// (it *averages* over orbit members); we rescale by `group.order()` so -/// that the result is the *summing* projector, matching -/// [`super::symmetry_merge_pauli_sum`]. +/// This is the **summing** projector +/// `Σ_{p ∈ orbit} χ_k(g_p) · c_p` over each orbit's *distinct* members, not the +/// orbit-averaged one that [`canonicalize_pauli_sum_complex`] returns. +/// Summing is what makes the merge idempotent — and hence safe to apply +/// after every Trotter step — for *every* orbit, including orbits with a +/// non-trivial stabilizer, and it reduces exactly to +/// [`super::symmetry_merge_pauli_sum`] at `k = 0`. /// -/// Entries whose rescaled component is exactly zero are dropped, so a -/// purely real operator leaves `im` empty. +/// Entries whose component is exactly zero are dropped, so a purely real +/// operator leaves `im` empty. /// /// # Panics /// @@ -208,24 +241,15 @@ pub fn momentum_merge_pauli_sum_pair( for (word, coeff) in im.data().iter() { combined.entry(*word).or_insert(Complex::new(0.0, 0.0)).im += *coeff; } - let mut basis = Vec::with_capacity(combined.len()); - let mut coeffs = Vec::with_capacity(combined.len()); - for (word, coeff) in combined { - basis.push(word); - coeffs.push(coeff); - } - - canonicalize_pauli_sum_complex(&mut basis, &mut coeffs, group, k_modes); - - let scale = group.order() as f64; + let projected = project_onto_reps(&combined, group, k_modes); re.data_mut().clear(); im.data_mut().clear(); - for (word, coeff) in basis.into_iter().zip(coeffs) { - if coeff.re != 0.0 { - *re += (word, coeff.re * scale); + for (word, (sum, _orbit_size)) in projected { + if sum.re != 0.0 { + *re += (word, sum.re); } - if coeff.im != 0.0 { - *im += (word, coeff.im * scale); + if sum.im != 0.0 { + *im += (word, sum.im); } } } diff --git a/crates/ppvm-pauli-sum/src/symmetry/tests.rs b/crates/ppvm-pauli-sum/src/symmetry/tests.rs index ee9981e60..8af9a01bf 100644 --- a/crates/ppvm-pauli-sum/src/symmetry/tests.rs +++ b/crates/ppvm-pauli-sum/src/symmetry/tests.rs @@ -473,10 +473,11 @@ where (re, im) } -/// At `k = 0` on free orbits, the phase-aware pair merge must agree -/// with the independent real-coefficient `symmetry_merge_pauli_sum` -/// code path — this is what pins the `|G|` rescale (the momentum -/// projector averages; `symmetry_merge` sums). +/// At `k = 0` the phase-aware pair merge must agree with the independent +/// real-coefficient `symmetry_merge_pauli_sum` code path — for *every* +/// orbit, free or stabilized. This is the +/// regression test for the summing-vs-averaging convention: a global +/// `|G|` rescale of the averaged projector agrees only on free orbits. #[test] fn momentum_merge_pair_matches_symmetry_merge_at_k0() { use crate::config::indexmap::ByteFxHashF64; @@ -486,10 +487,10 @@ fn momentum_merge_pair_matches_symmetry_merge_at_k0() { let n = 4usize; let group = TranslationGroup::chain_1d(n); - // Σ_j Z_j and Σ_j X_j X_{j+1}: both free orbits (|orbit| = |G|). let mut reference: PauliSum = PauliSum::builder().n_qubits(n).build(); let mut re: PauliSum = PauliSum::builder().n_qubits(n).build(); let mut im: PauliSum = PauliSum::builder().n_qubits(n).build(); + // Σ_j Z_j and Σ_j X_j X_{j+1}: free orbits (|orbit| = |G|). for j in 0..n { let mut z: Vec = vec!['I'; n]; z[j] = 'Z'; @@ -502,6 +503,12 @@ fn momentum_merge_pair_matches_symmetry_merge_at_k0() { re += (st.as_str(), c); } } + // Orbits WITH a stabilizer, where a global |G| rescale would be wrong: + // "ZZZZ" is translation-invariant (|orbit| = 1) and "ZIZI" has period 2. + for (st, c) in [("ZZZZ", 0.75), ("ZIZI", 2.0), ("IZIZ", -0.5)] { + reference += (st, c); + re += (st, c); + } // `im` needs an entry to exist; a zero coefficient must not survive. im += ("IIII", 0.0); diff --git a/crates/ppvm-python-native/src/interface.rs b/crates/ppvm-python-native/src/interface.rs index e84c29fe5..7949b9b24 100644 --- a/crates/ppvm-python-native/src/interface.rs +++ b/crates/ppvm-python-native/src/interface.rs @@ -91,8 +91,9 @@ macro_rules! create_interface_symmetry_methods { /// reduces to `symmetry_merge`). This generalizes /// `symmetry_merge` to k != 0 while keeping real coefficients on /// the Python side — the only place complex arithmetic appears - /// is the internal character-weighted fold, reusing the tested - /// `canonicalize_pauli_sum_complex`. + /// is the internal character-weighted fold, which like + /// `symmetry_merge` *sums* over each orbit (so the merge is + /// idempotent on every orbit, free or stabilized). /// /// `self` and `other` must be distinct objects with identical /// qubit count. After a translation-covariant gate layer this diff --git a/crates/ppvm-python-native/src/symmetry.rs b/crates/ppvm-python-native/src/symmetry.rs index 880851468..bd03a1b72 100644 --- a/crates/ppvm-python-native/src/symmetry.rs +++ b/crates/ppvm-python-native/src/symmetry.rs @@ -171,8 +171,9 @@ impl TranslationGroup { /// `momentum` is a length-`group.n_generators` integer array of mode /// indices; the wavenumber along generator `g` is /// `2π · momentum[g] / group.generator_order(g)`. Use `momentum=[0, …]` -/// for the trivial (k=0) sector — equivalent to plain merging modulo -/// the 1/|G| normalization the complex merge applies. +/// for the trivial (k=0) sector — equivalent to plain merging modulo the +/// `1/|orbit|` normalization this projection applies (it *averages* over +/// each orbit's distinct members; `PauliSum.momentum_merge` sums). /// /// If the input is **not** in sector `momentum`, the projection /// silently throws away the other components. Use diff --git a/ppvm-python/src/ppvm/paulisum.py b/ppvm-python/src/ppvm/paulisum.py index b367a7244..3504cd712 100644 --- a/ppvm-python/src/ppvm/paulisum.py +++ b/ppvm-python/src/ppvm/paulisum.py @@ -416,7 +416,12 @@ def momentum_merge(self, other: "PauliSum", group, momentum) -> None: Generalizes `symmetry_merge` to non-trivial momentum sectors (``k != 0``) while keeping real coefficients on both PauliSums — the - only complex arithmetic is the internal character-weighted fold. + only complex arithmetic is the internal character-weighted fold. Like + `symmetry_merge` it is the *summing* projector + ``Σ_{p in orbit} χ_k(g_p)·c_p``, hence idempotent on every orbit and + safe to apply after each Trotter step; at ``momentum=[0, ...]`` it + reduces exactly to `symmetry_merge`. + ``self`` and ``other`` must be distinct objects with the same qubit count. Exact after a translation-covariant gate layer; under a generic Trotter step it carries the same ``O(dt^{p+1})`` equivariance diff --git a/ppvm-python/test/test_momentum_merge.py b/ppvm-python/test/test_momentum_merge.py index 6dccfc84d..8e5fb953f 100644 --- a/ppvm-python/test/test_momentum_merge.py +++ b/ppvm-python/test/test_momentum_merge.py @@ -91,19 +91,16 @@ def test_momentum_merge_idempotent(k): assert max(abs(once.get(x, 0j) - twice.get(x, 0j)) for x in keys) < 1e-12 -@pytest.mark.xfail( - reason="momentum_merge rescales by |G| but the projector averages over the " - "|orbit| DISTINCT members, so orbits with a non-trivial stabilizer are " - "amplified by |G|/|orbit| on every merge. `_seed_pair` only produces free " - "orbits, which is why test_momentum_merge_idempotent passes.", - strict=True, -) -@pytest.mark.parametrize("word, orbit_size", [("ZZZZ", 1), ("ZIZI", 2)]) -def test_momentum_merge_idempotent_on_stabilized_orbit(word, orbit_size): +@pytest.mark.parametrize("word", ["ZZZZ", "ZIZI"]) +def test_momentum_merge_idempotent_on_stabilized_orbit(word): """Idempotency must hold for every orbit, not just the free ones. - ``ZZZZ`` is translation-invariant (orbit size 1) and ``ZIZI`` has period 2, - so on a 4-site chain they pick up factors of 4 and 2 per merge. + ``ZZZZ`` is translation-invariant (orbit size 1) and ``ZIZI`` has period + 2, so on a 4-site chain their orbits are smaller than the group. A merge + that rescaled the orbit-*averaged* projection by a global ``|G|`` would + amplify them by ``|G|/|orbit|`` — 4x and 2x — on every merge; the summing + projector leaves them fixed. ``_seed_pair`` only produces free orbits, + where the two conventions coincide, so this case needs its own test. """ n = 4 g = TranslationGroup.chain_1d(n) @@ -111,6 +108,8 @@ def test_momentum_merge_idempotent_on_stabilized_orbit(word, orbit_size): PB = PauliSum.new(n, [(word, 0.0)], min_abs_coeff=0.0, max_pauli_weight=n) PA.momentum_merge(PB, g, [0]) once = _to_complex_dict(PA, PB) + # The coefficient is conserved outright, not just stable under re-merging. + assert sorted(abs(v) for v in once.values() if abs(v) > 1e-12) == pytest.approx([1.0]) PA.momentum_merge(PB, g, [0]) twice = _to_complex_dict(PA, PB) keys = set(once) | set(twice) From 74327041901854f348a2c8ce4e609147e6af34b2 Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 2 Sep 2026 14:14:27 +0200 Subject: [PATCH 11/13] refactor(python-native): share the Pauli-array codec; wrap the symmetry API MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two related pieces of boilerplate cleanup on the `(basis_arr, coeffs)` array surface. Rust: the validate-decode-encode triple was written out three times in `symmetry.rs` and again in `lindblad.rs::pc_step_orbit_rep`, and `symmetry.rs` reached across into `crate::lindblad::decode_basis` for the decode half. All of it now lives in one `pauli_arr` module — `decode_basis` (moved), `encode_basis`, and the three argument checks (`check_group_width`, `check_coeffs_len`, `check_momentum_len`). Error messages are unchanged, so the existing tests that match on them still pass. Net -155 lines across `lindblad.rs` and `symmetry.rs`. One incidental improvement: `check_momentum_sector_arr` never validated its `coeffs` / `momentum` lengths, so a mismatch reached the core's `assert_eq!` and surfaced as a PanicException. It now raises ValueError like its siblings. Python: the three `_core` symmetry functions had no wrapper, so they demanded exact dtypes and rejected the natural `np.array([1])` momentum (numpy's default int64) with `TypeError: 'ndarray' object is not an instance of 'ndarray'`. A new `ppvm.symmetry` module wraps them with the same dtype coercion `lindblad.py` already applies, and re-exports `TranslationGroup` — which previously had to be imported from the private `ppvm._core` despite the docstrings pointing users at it. `ppvm` now exports `TranslationGroup`, `canonicalize_basis_arr`, `canonicalize_basis_arr_complex` and `check_momentum_sector_arr`; the docstrings that said `ppvm._core.TranslationGroup` are updated, and griffe picks the module up for the docs site automatically. 285 Python tests pass (one new, for the coercion), `cargo test --workspace` is green, clippy/ruff/ty are clean. Co-Authored-By: Claude Opus 5 (1M context) --- crates/ppvm-python-native/src/interface.rs | 2 +- crates/ppvm-python-native/src/lib.rs | 1 + crates/ppvm-python-native/src/lindblad.rs | 71 +------- crates/ppvm-python-native/src/pauli_arr.rs | 94 +++++++++++ crates/ppvm-python-native/src/symmetry.rs | 84 ++-------- ppvm-python/src/ppvm/__init__.py | 4 + ppvm-python/src/ppvm/paulisum.py | 4 +- ppvm-python/src/ppvm/symmetry.py | 158 ++++++++++++++++++ .../test/lindblad/test_pc_step_orbit_rep.py | 7 +- ppvm-python/test/test_momentum_merge.py | 3 +- ppvm-python/test/test_symmetry_arrays.py | 34 +++- ppvm-python/test/test_symmetry_merge.py | 3 +- 12 files changed, 322 insertions(+), 143 deletions(-) create mode 100644 crates/ppvm-python-native/src/pauli_arr.rs create mode 100644 ppvm-python/src/ppvm/symmetry.py diff --git a/crates/ppvm-python-native/src/interface.rs b/crates/ppvm-python-native/src/interface.rs index 7949b9b24..20450252b 100644 --- a/crates/ppvm-python-native/src/interface.rs +++ b/crates/ppvm-python-native/src/interface.rs @@ -62,7 +62,7 @@ macro_rules! create_interface_symmetry_methods { /// entry count by up to `|group|×` for translation-invariant /// operators. /// - /// See `ppvm._core.TranslationGroup` for constructors + /// See `ppvm.TranslationGroup` for constructors /// (`chain_1d`, `torus_2d`, `torus_3d`, `ladder`). /// /// Plain real-coefficient merge (the `k=0` symmetry sector). diff --git a/crates/ppvm-python-native/src/lib.rs b/crates/ppvm-python-native/src/lib.rs index efb190104..e4f9228c8 100644 --- a/crates/ppvm-python-native/src/lib.rs +++ b/crates/ppvm-python-native/src/lib.rs @@ -15,6 +15,7 @@ pub mod interface; pub mod interface_tableau; pub mod interface_tableau_sum; pub mod lindblad; +pub mod pauli_arr; pub mod stim_program; pub mod symmetry; diff --git a/crates/ppvm-python-native/src/lindblad.rs b/crates/ppvm-python-native/src/lindblad.rs index 88799a321..e56b78351 100644 --- a/crates/ppvm-python-native/src/lindblad.rs +++ b/crates/ppvm-python-native/src/lindblad.rs @@ -13,10 +13,8 @@ use std::collections::HashMap; use num::Complex; -use numpy::{ - Complex64, IntoPyArray, PyArray1, PyArray2, PyArrayMethods, PyReadonlyArray1, PyReadonlyArray2, -}; -use ppvm_lindblad::{JumpInput, LindbladSpec as CoreSpec, Word, codes_from_word, word_from_codes}; +use numpy::{Complex64, IntoPyArray, PyArray1, PyArray2, PyReadonlyArray1, PyReadonlyArray2}; +use ppvm_lindblad::{JumpInput, LindbladSpec as CoreSpec, Word, word_from_codes}; use pyo3::{exceptions::PyValueError, prelude::*}; type PyPauliMap<'py> = (Bound<'py, PyArray2>, Bound<'py, PyArray1>); @@ -27,7 +25,7 @@ type PyCoo<'py> = ( Bound<'py, PyArray1>, ); -fn map_err(e: ppvm_lindblad::Error) -> PyErr { +pub(crate) fn map_err(e: ppvm_lindblad::Error) -> PyErr { PyValueError::new_err(e.to_string()) } @@ -46,29 +44,7 @@ fn assert_basis_unique(basis: &[Word]) -> PyResult<()> { Ok(()) } -/// Decode a `(N, n_qubits)` uint8 ndarray view into `N` packed [`Word`]s. -pub(crate) fn decode_basis( - view: &numpy::ndarray::ArrayView2, - n_qubits: usize, -) -> PyResult> { - let n_basis = view.shape()[0]; - let n_cols = view.shape()[1]; - if n_cols != n_qubits { - return Err(PyValueError::new_err(format!( - "basis has {n_cols} columns but spec.n_qubits = {n_qubits}" - ))); - } - let mut out = Vec::with_capacity(n_basis); - let mut row_buf = vec![0u8; n_qubits]; - for i in 0..n_basis { - let row = view.row(i); - for (q, slot) in row_buf.iter_mut().enumerate() { - *slot = row[q]; - } - out.push(word_from_codes(&row_buf).map_err(map_err)?); - } - Ok(out) -} +use crate::pauli_arr::{check_coeffs_len, check_momentum_len, decode_basis, encode_basis}; /// Pack `Vec<(Word, f64)>` into the standard PyO3 return shape. fn pack_pauli_map<'py>( @@ -76,17 +52,8 @@ fn pack_pauli_map<'py>( pairs: Vec<(Word, f64)>, n_qubits: usize, ) -> PyResult> { - let m = pairs.len(); - let mut basis = vec![0u8; m * n_qubits]; - let mut coeffs = vec![0f64; m]; - for (i, (w, c)) in pairs.into_iter().enumerate() { - codes_from_word(&w, &mut basis[i * n_qubits..(i + 1) * n_qubits]); - coeffs[i] = c; - } - let basis_arr = basis - .into_pyarray(py) - .reshape([m, n_qubits]) - .map_err(|e| PyValueError::new_err(format!("reshape failed: {e}")))?; + let (words, coeffs): (Vec, Vec) = pairs.into_iter().unzip(); + let basis_arr = encode_basis(py, &words, n_qubits)?; Ok((basis_arr, coeffs.into_pyarray(py))) } @@ -178,13 +145,7 @@ impl LindbladSpec { let basis_view = basis.as_array(); let basis_words = decode_basis(&basis_view, n_q)?; let coeffs_slice = coeffs.as_slice()?; - if coeffs_slice.len() != basis_words.len() { - return Err(PyValueError::new_err(format!( - "coeffs has length {} but basis has {} rows", - coeffs_slice.len(), - basis_words.len() - ))); - } + check_coeffs_len(coeffs_slice.len(), basis_words.len())?; let protected_words: Vec = if let Some(ref prot) = protected { let pv = prot.as_array(); decode_basis(&pv, n_q)? @@ -414,13 +375,7 @@ impl LindbladSpec { Vec::new() }; let k_slice = momentum.as_slice()?; - if k_slice.len() != group.core().n_generators() { - return Err(PyValueError::new_err(format!( - "momentum has {} entries but group has {} generators", - k_slice.len(), - group.core().n_generators() - ))); - } + check_momentum_len(k_slice.len(), group.core().n_generators())?; if canonicalize_first { canonicalize_basis_to_rep(&mut basis_words, group.core()); } @@ -441,19 +396,11 @@ impl LindbladSpec { ) .map_err(map_err)?; - let m = basis_words.len(); - let mut out_basis = vec![0u8; m * n_q]; - for (i, w) in basis_words.iter().enumerate() { - codes_from_word(w, &mut out_basis[i * n_q..(i + 1) * n_q]); - } let out_coeffs: Vec = coeffs_vec .iter() .map(|c| Complex64::new(c.re, c.im)) .collect(); - let basis_arr = out_basis - .into_pyarray(py) - .reshape([m, n_q]) - .map_err(|e| PyValueError::new_err(format!("reshape failed: {e}")))?; + let basis_arr = encode_basis(py, &basis_words, n_q)?; Ok((basis_arr, out_coeffs.into_pyarray(py))) } diff --git a/crates/ppvm-python-native/src/pauli_arr.rs b/crates/ppvm-python-native/src/pauli_arr.rs new file mode 100644 index 000000000..a90402cc4 --- /dev/null +++ b/crates/ppvm-python-native/src/pauli_arr.rs @@ -0,0 +1,94 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Shared codec and argument validation for the `(N, n_qubits)` uint8 +//! Pauli-basis array representation used across the Lindblad and +//! symmetry bindings. +//! +//! Every `*_arr` entry point decodes an incoming basis array into packed +//! [`Word`]s, validates the companion `coeffs` / `momentum` lengths, and +//! re-encodes the result on the way out. These are those three steps. + +use numpy::{IntoPyArray, PyArray2, PyArrayMethods}; +use ppvm_lindblad::{Word, codes_from_word, word_from_codes}; +use pyo3::{exceptions::PyValueError, prelude::*}; + +use crate::lindblad::map_err; + +/// Decode a `(N, n_qubits)` uint8 ndarray view into `N` packed [`Word`]s. +pub(crate) fn decode_basis( + view: &numpy::ndarray::ArrayView2, + n_qubits: usize, +) -> PyResult> { + let n_basis = view.shape()[0]; + let n_cols = view.shape()[1]; + if n_cols != n_qubits { + return Err(PyValueError::new_err(format!( + "basis has {n_cols} columns but spec.n_qubits = {n_qubits}" + ))); + } + let mut out = Vec::with_capacity(n_basis); + let mut row_buf = vec![0u8; n_qubits]; + for i in 0..n_basis { + let row = view.row(i); + for (q, slot) in row_buf.iter_mut().enumerate() { + *slot = row[q]; + } + out.push(word_from_codes(&row_buf).map_err(map_err)?); + } + Ok(out) +} + +/// Encode packed [`Word`]s back into an `(M, n_qubits)` uint8 array. +pub(crate) fn encode_basis<'py>( + py: Python<'py>, + words: &[Word], + n_qubits: usize, +) -> PyResult>> { + let m = words.len(); + let mut flat = vec![0u8; m * n_qubits]; + for (i, w) in words.iter().enumerate() { + codes_from_word(w, &mut flat[i * n_qubits..(i + 1) * n_qubits]); + } + flat.into_pyarray(py) + .reshape([m, n_qubits]) + .map_err(|e| PyValueError::new_err(format!("reshape failed: {e}"))) +} + +/// Check the row width of a basis array against the qubit count a +/// [`crate::symmetry::TranslationGroup`] acts on. Reported separately from +/// [`decode_basis`]'s own width check so the error names the group rather +/// than the spec. +pub(crate) fn check_group_width( + view: &numpy::ndarray::ArrayView2, + n_qubits: usize, +) -> PyResult<()> { + let width = view.shape().get(1).copied(); + if width != Some(n_qubits) { + return Err(PyValueError::new_err(format!( + "basis has {} qubits per row but group acts on {n_qubits}", + width.unwrap_or(0) + ))); + } + Ok(()) +} + +/// Check that a coefficient vector has one entry per basis row. +pub(crate) fn check_coeffs_len(n_coeffs: usize, n_rows: usize) -> PyResult<()> { + if n_coeffs != n_rows { + return Err(PyValueError::new_err(format!( + "coeffs has length {n_coeffs} but basis has {n_rows} rows" + ))); + } + Ok(()) +} + +/// Check that a momentum vector has one mode index per group generator. +pub(crate) fn check_momentum_len(n_modes: usize, n_generators: usize) -> PyResult<()> { + if n_modes != n_generators { + return Err(PyValueError::new_err(format!( + "momentum has {n_modes} entries but group has {n_generators} generators" + ))); + } + Ok(()) +} diff --git a/crates/ppvm-python-native/src/symmetry.rs b/crates/ppvm-python-native/src/symmetry.rs index bd03a1b72..6c060ec51 100644 --- a/crates/ppvm-python-native/src/symmetry.rs +++ b/crates/ppvm-python-native/src/symmetry.rs @@ -11,13 +11,15 @@ //! used by `Lindbladian.pc_step_arr`. use num::Complex; -use numpy::{ - Complex64, IntoPyArray, PyArray1, PyArray2, PyArrayMethods, PyReadonlyArray1, PyReadonlyArray2, -}; +use numpy::{Complex64, IntoPyArray, PyArray1, PyArray2, PyReadonlyArray1, PyReadonlyArray2}; use ppvm_lindblad::{codes_from_word, word_from_codes}; use ppvm_pauli_sum::symmetry as core_sym; use pyo3::{exceptions::PyValueError, prelude::*}; +use crate::pauli_arr::{ + check_coeffs_len, check_group_width, check_momentum_len, decode_basis, encode_basis, +}; + type PyPauliMap<'py> = (Bound<'py, PyArray2>, Bound<'py, PyArray1>); type PyPauliMapComplex<'py> = (Bound<'py, PyArray2>, Bound<'py, PyArray1>); @@ -188,30 +190,12 @@ pub fn canonicalize_basis_arr_complex<'py>( ) -> PyResult> { let basis_view = basis.as_array(); let n_q = group.inner.n_qubits(); - if basis_view.shape().get(1).copied() != Some(n_q) { - return Err(PyValueError::new_err(format!( - "basis has {} qubits per row but group acts on {n_q}", - basis_view.shape().get(1).copied().unwrap_or(0) - ))); - } - let n = basis_view.shape()[0]; + check_group_width(&basis_view, n_q)?; let coeffs_slice = coeffs.as_slice()?; - if coeffs_slice.len() != n { - return Err(PyValueError::new_err(format!( - "coeffs has length {} but basis has {} rows", - coeffs_slice.len(), - n - ))); - } + check_coeffs_len(coeffs_slice.len(), basis_view.shape()[0])?; let k_slice = momentum.as_slice()?; - if k_slice.len() != group.inner.n_generators() { - return Err(PyValueError::new_err(format!( - "momentum has {} entries but group has {} generators", - k_slice.len(), - group.inner.n_generators() - ))); - } - let mut basis_words = crate::lindblad::decode_basis(&basis_view, n_q)?; + check_momentum_len(k_slice.len(), group.inner.n_generators())?; + let mut basis_words = decode_basis(&basis_view, n_q)?; let mut coeffs_vec: Vec> = coeffs_slice .iter() .map(|c| Complex::new(c.re, c.im)) @@ -224,19 +208,11 @@ pub fn canonicalize_basis_arr_complex<'py>( k_slice, ); - let m = basis_words.len(); - let mut out_basis = vec![0u8; m * n_q]; - for (i, w) in basis_words.iter().enumerate() { - codes_from_word(w, &mut out_basis[i * n_q..(i + 1) * n_q]); - } let out_coeffs: Vec = coeffs_vec .iter() .map(|c| Complex64::new(c.re, c.im)) .collect(); - let basis_arr = out_basis - .into_pyarray(py) - .reshape([m, n_q]) - .map_err(|e| PyValueError::new_err(format!("reshape failed: {e}")))?; + let basis_arr = encode_basis(py, &basis_words, n_q)?; Ok((basis_arr, out_coeffs.into_pyarray(py))) } @@ -257,15 +233,12 @@ pub fn check_momentum_sector_arr<'py>( ) -> PyResult<()> { let basis_view = basis.as_array(); let n_q = group.inner.n_qubits(); - if basis_view.shape().get(1).copied() != Some(n_q) { - return Err(PyValueError::new_err(format!( - "basis has {} qubits per row but group acts on {n_q}", - basis_view.shape().get(1).copied().unwrap_or(0) - ))); - } + check_group_width(&basis_view, n_q)?; let coeffs_slice = coeffs.as_slice()?; + check_coeffs_len(coeffs_slice.len(), basis_view.shape()[0])?; let k_slice = momentum.as_slice()?; - let basis_words = crate::lindblad::decode_basis(&basis_view, n_q)?; + check_momentum_len(k_slice.len(), group.inner.n_generators())?; + let basis_words = decode_basis(&basis_view, n_q)?; let coeffs_vec: Vec> = coeffs_slice .iter() .map(|c| Complex::new(c.re, c.im)) @@ -295,36 +268,15 @@ pub fn canonicalize_basis_arr<'py>( ) -> PyResult> { let basis_view = basis.as_array(); let n_q = group.inner.n_qubits(); - if basis_view.shape().get(1).copied() != Some(n_q) { - return Err(PyValueError::new_err(format!( - "basis has {} qubits per row but group acts on {n_q}", - basis_view.shape().get(1).copied().unwrap_or(0) - ))); - } - let n = basis_view.shape()[0]; + check_group_width(&basis_view, n_q)?; let coeffs_slice = coeffs.as_slice()?; - if coeffs_slice.len() != n { - return Err(PyValueError::new_err(format!( - "coeffs has length {} but basis has {} rows", - coeffs_slice.len(), - n - ))); - } + check_coeffs_len(coeffs_slice.len(), basis_view.shape()[0])?; - let mut basis_words = crate::lindblad::decode_basis(&basis_view, n_q)?; + let mut basis_words = decode_basis(&basis_view, n_q)?; let mut coeffs_vec = coeffs_slice.to_vec(); core_sym::canonicalize_pauli_sum(&mut basis_words, &mut coeffs_vec, &group.inner); - // Re-encode. - let m = basis_words.len(); - let mut out_basis = vec![0u8; m * n_q]; - for (i, w) in basis_words.iter().enumerate() { - codes_from_word(w, &mut out_basis[i * n_q..(i + 1) * n_q]); - } - let basis_arr = out_basis - .into_pyarray(py) - .reshape([m, n_q]) - .map_err(|e| PyValueError::new_err(format!("reshape failed: {e}")))?; + let basis_arr = encode_basis(py, &basis_words, n_q)?; Ok((basis_arr, coeffs_vec.into_pyarray(py))) } diff --git a/ppvm-python/src/ppvm/__init__.py b/ppvm-python/src/ppvm/__init__.py index 286250622..7c1b0fa8b 100644 --- a/ppvm-python/src/ppvm/__init__.py +++ b/ppvm-python/src/ppvm/__init__.py @@ -16,3 +16,7 @@ from .squin_interpreter.device import ( GeneralizedTableauSimulatorTask as GeneralizedTableauSimulatorTask, ) +from .symmetry import TranslationGroup as TranslationGroup +from .symmetry import canonicalize_basis_arr as canonicalize_basis_arr +from .symmetry import canonicalize_basis_arr_complex as canonicalize_basis_arr_complex +from .symmetry import check_momentum_sector_arr as check_momentum_sector_arr diff --git a/ppvm-python/src/ppvm/paulisum.py b/ppvm-python/src/ppvm/paulisum.py index 3504cd712..f71c6d846 100644 --- a/ppvm-python/src/ppvm/paulisum.py +++ b/ppvm-python/src/ppvm/paulisum.py @@ -401,7 +401,7 @@ def symmetry_merge(self, group) -> None: sector only. Args: - group: A `ppvm._core.TranslationGroup` + group: A `ppvm.TranslationGroup` (use ``TranslationGroup.chain_1d(n)``, ``.torus_2d``, ``.torus_3d``, ``.ladder``, or ``.from_generators``). """ @@ -429,7 +429,7 @@ def momentum_merge(self, other: "PauliSum", group, momentum) -> None: Args: other: the PauliSum holding the imaginary component (modified in place). - group: a `ppvm._core.TranslationGroup`. + group: a `ppvm.TranslationGroup`. momentum: sequence of integer modes, one per group generator (e.g. ``[k]`` for a 1D chain; ``[0, ...]`` is the trivial sector). """ diff --git a/ppvm-python/src/ppvm/symmetry.py b/ppvm-python/src/ppvm/symmetry.py new file mode 100644 index 000000000..d7a7ee9b4 --- /dev/null +++ b/ppvm-python/src/ppvm/symmetry.py @@ -0,0 +1,158 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 + +"""Translation-symmetry merging of Pauli sums. + +A `TranslationGroup` is a finite abelian group acting on qubit positions by +permutation. Every Pauli word then belongs to a translation orbit, and +dynamics that commutes with the group can be tracked using **one canonical +representative per orbit** instead of all ``|G|`` members — cutting per-step +memory and compute by up to ``|G|×`` (Teng et al., arXiv:2512.12094). + +Two representations are supported: + +- `ppvm.PauliSum.symmetry_merge` / `ppvm.PauliSum.momentum_merge` for the + dictionary representation used by gate-based Trotter evolution. +- the ``*_arr`` functions here for the ``(basis_arr, coeffs)`` array + representation used by `ppvm.Lindbladian.pc_step_arr` and + `ppvm.Lindbladian.pc_step_orbit_rep`, where ``basis_arr`` is an + ``(N, n_qubits)`` uint8 array with the encoding ``0=I, 1=X, 2=Z, 3=Y``. + +These are thin wrappers over `ppvm._core` that coerce their arguments to the +dtypes the compiled entry points require (uint8 basis, float64 / complex128 +coefficients, int32 momentum), so plain Python lists and default-dtype numpy +arrays work. +""" + +from __future__ import annotations + +import numpy as np +import numpy.typing as npt + +from . import _core +from ._core import TranslationGroup as TranslationGroup + +__all__ = [ + "TranslationGroup", + "canonicalize_basis_arr", + "canonicalize_basis_arr_complex", + "check_momentum_sector_arr", +] + + +def _momentum(momentum: npt.ArrayLike) -> np.ndarray: + return np.ascontiguousarray(momentum, dtype=np.int32) + + +def _basis(basis_arr: npt.ArrayLike) -> np.ndarray: + return np.ascontiguousarray(basis_arr, dtype=np.uint8) + + +def canonicalize_basis_arr( + basis_arr: npt.ArrayLike, + coeffs: npt.ArrayLike, + group: TranslationGroup, +) -> tuple[np.ndarray, np.ndarray]: + """Merge a real-coefficient ``(basis_arr, coeffs)`` Pauli sum into + orbit-representative form. + + Each row of ``basis_arr`` is replaced by its canonical representative + under ``group``; coefficients of rows collapsing to the same + representative are **summed**. The output is no longer than the input. + + This is the trivial (``k=0``) symmetry sector. For dynamics that commutes + with ``group`` and a ``group``-invariant initial state, it preserves every + ``group``-invariant expectation value. Use `canonicalize_basis_arr_complex` + for non-trivial momentum sectors. + + Args: + basis_arr: ``(N, n_qubits)`` array of Pauli codes. + coeffs: length-``N`` real coefficients. + group: the symmetry group to merge under. + + Returns: + ``(merged_basis_arr, merged_coeffs)``. + """ + return _core.canonicalize_basis_arr( + _basis(basis_arr), + np.ascontiguousarray(coeffs, dtype=np.float64), + group, + ) + + +def canonicalize_basis_arr_complex( + basis_arr: npt.ArrayLike, + coeffs: npt.ArrayLike, + group: TranslationGroup, + momentum: npt.ArrayLike, +) -> tuple[np.ndarray, np.ndarray]: + """Phase-aware merge of a complex-coefficient ``(basis_arr, coeffs)`` + Pauli sum into orbit-representative form, projected onto momentum sector + ``momentum``. + + Coefficients on each orbit's distinct members are **averaged** with the + character weight ``χ_k(g)`` — a ``1/|orbit|`` normalization that + `canonicalize_basis_arr` (which sums) does not apply. Orbits whose + stabilizer cannot carry ``momentum`` project to zero and are dropped. + + If the input does not actually lie in sector ``momentum``, the projection + silently discards the other components; call `check_momentum_sector_arr` + first to validate. + + Args: + basis_arr: ``(N, n_qubits)`` array of Pauli codes. + coeffs: length-``N`` complex coefficients. + group: the symmetry group to merge under. + momentum: one integer mode index per group generator. The wavenumber + along generator ``g`` is ``2π · momentum[g] / order_g``; + ``[0, ...]`` is the trivial sector. + + Returns: + ``(merged_basis_arr, merged_coeffs)`` with complex coefficients. + """ + return _core.canonicalize_basis_arr_complex( + _basis(basis_arr), + np.ascontiguousarray(coeffs, dtype=np.complex128), + group, + _momentum(momentum), + ) + + +def check_momentum_sector_arr( + basis_arr: npt.ArrayLike, + coeffs: npt.ArrayLike, + group: TranslationGroup, + momentum: npt.ArrayLike, + tol: float = 1e-8, +) -> None: + """Verify that a ``(basis_arr, complex_coeffs)`` Pauli sum lies entirely + in momentum sector ``momentum``. + + For every orbit represented in the basis, all members must satisfy + ``c_{g·r} = χ_k(g)⁻¹ · c_r``. Orbit members absent from ``basis_arr`` + count as zero rather than being ignored, so a partially-populated orbit + fails. + + Run this on a user-supplied initial state before feeding it to + `canonicalize_basis_arr_complex` or + `ppvm.Lindbladian.pc_step_orbit_rep` — silently projecting a + wrongly-typed input throws away meaningful physics. + + Args: + basis_arr: ``(N, n_qubits)`` array of Pauli codes. + coeffs: length-``N`` complex coefficients. + group: the symmetry group. + momentum: one integer mode index per group generator. + tol: relative tolerance on the coefficient comparison. + + Raises: + ValueError: if the input is not in the sector, naming the offending + orbit representative with its expected and actual coefficient. + """ + return _core.check_momentum_sector_arr( + _basis(basis_arr), + np.ascontiguousarray(coeffs, dtype=np.complex128), + group, + _momentum(momentum), + tol, + ) diff --git a/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py b/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py index fdb234cdb..36dddaf1b 100644 --- a/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py +++ b/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py @@ -20,8 +20,7 @@ import numpy as np import pytest -from ppvm import Lindbladian -from ppvm._core import TranslationGroup, canonicalize_basis_arr_complex +from ppvm import Lindbladian, TranslationGroup, canonicalize_basis_arr_complex from ._helpers import all_strings @@ -45,7 +44,9 @@ def to_dict(basis, coeffs): def momentum(*modes): - return np.array(modes, dtype=np.int32) + """A plain tuple: both the wrappers and `Lindbladian.pc_step_orbit_rep` + coerce momentum to the int32 the compiled code needs.""" + return modes def xy_chain_pbc(n, gamma): diff --git a/ppvm-python/test/test_momentum_merge.py b/ppvm-python/test/test_momentum_merge.py index 8e5fb953f..217863299 100644 --- a/ppvm-python/test/test_momentum_merge.py +++ b/ppvm-python/test/test_momentum_merge.py @@ -18,8 +18,7 @@ import numpy as np import pytest -from ppvm import PauliSum -from ppvm._core import TranslationGroup +from ppvm import PauliSum, TranslationGroup # ── dense Pauli helpers (exact references) ─────────────────────────────────── _I = np.eye(2, dtype=complex) diff --git a/ppvm-python/test/test_symmetry_arrays.py b/ppvm-python/test/test_symmetry_arrays.py index 11ec2aa59..7630e6cc0 100644 --- a/ppvm-python/test/test_symmetry_arrays.py +++ b/ppvm-python/test/test_symmetry_arrays.py @@ -2,7 +2,8 @@ # SPDX-License-Identifier: Apache-2.0 """Tests for the array-form symmetry primitives on the ``(basis_arr, coeffs)`` -representation used by ``Lindbladian.pc_step_arr``: +representation used by ``Lindbladian.pc_step_arr``, as exported from the +``ppvm`` package (thin dtype-coercing wrappers over ``ppvm._core``): - ``canonicalize_basis_arr`` — plain real merge (sums colliding coefficients) - ``canonicalize_basis_arr_complex`` — momentum-sector projection (averages @@ -19,7 +20,7 @@ import numpy as np import pytest -from ppvm._core import ( +from ppvm import ( TranslationGroup, canonicalize_basis_arr, canonicalize_basis_arr_complex, @@ -48,15 +49,38 @@ def rep_of(group, s): def momentum(*modes): - """The ``_core`` free functions take momentum as an int32 array; numpy's - default integer dtype is int64, which they reject.""" - return np.array(modes, dtype=np.int32) + """The wrappers coerce momentum for us; most tests pass a plain tuple. + + See `test_wrappers_coerce_argument_dtypes` for the coercion itself. + """ + return modes def z_strings(n): return ["I" * j + "Z" + "I" * (n - j - 1) for j in range(n)] +# ── argument coercion (the reason the wrappers exist) ──────────────────────── +def test_wrappers_coerce_argument_dtypes(): + """The compiled entry points demand exact dtypes — uint8 basis, float64 / + complex128 coefficients, int32 momentum. numpy's default integer dtype is + int64, so an unwrapped ``np.array([0])`` momentum is rejected; the wrappers + accept plain Python sequences and default-dtype arrays. + """ + n = 4 + g = TranslationGroup.chain_1d(n) + words = z_strings(n) + py_basis = [[_CODE[c] for c in s] for s in words] # list[list[int]] + + real = to_dict(canonicalize_basis_arr(py_basis, [1.0] * n, g)) + assert real == pytest.approx({rep_of(g, words[0]): float(n)}) + + # int64 momentum (numpy default) and a plain list of complex. + cx = to_dict(canonicalize_basis_arr_complex(py_basis, [1 + 0j] * n, g, np.array([0]))) + assert len(cx) == 1 + assert check_momentum_sector_arr(py_basis, [1 + 0j] * n, g, [0]) is None + + # ── canonicalize_basis_arr (real, k=0) ─────────────────────────────────────── def test_canonicalize_basis_arr_sums_collisions(): n = 4 diff --git a/ppvm-python/test/test_symmetry_merge.py b/ppvm-python/test/test_symmetry_merge.py index d151cc168..4de560263 100644 --- a/ppvm-python/test/test_symmetry_merge.py +++ b/ppvm-python/test/test_symmetry_merge.py @@ -12,8 +12,7 @@ import numpy as np import pytest -from ppvm import PauliSum -from ppvm._core import TranslationGroup +from ppvm import PauliSum, TranslationGroup _CODE = {"I": 0, "X": 1, "Z": 2, "Y": 3} _CHAR = {v: k for k, v in _CODE.items()} From 3a56041c03e037ef6c36d47875a80fbe745e7d9e Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 2 Sep 2026 14:23:01 +0200 Subject: [PATCH 12/13] fix(lindblad): honour num_threads on the orbit-rep step; API polish MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three small fixes closing out the symmetry-evolution cleanup. `pc_step_orbit_rep` silently ignored `cfg.num_threads` while its sibling `pc_step` honoured it — the orbit-rep path never went through `run_in_pool`. It does now, and `num_threads` is threaded through the PyO3 signature, the `.pyi` stub and the `lindblad.py` wrapper so the two adaptive step entry points take the same knobs. Previously the binding hardcoded `num_threads: None`, so the argument could not even be passed. `PauliSum.momentum_merge` took `other` as `PyRefMut`, so passing the same object as both the real and imaginary part died with PyO3's raw already-borrowed error instead of the "must be distinct objects" the docstring promised. It now takes a `Bound` and converts the failed borrow into that message. `group` / `momentum` parameters on `PauliSum.symmetry_merge`, `PauliSum.momentum_merge` and `Lindbladian.pc_step_orbit_rep` were unannotated; they now carry `_core.TranslationGroup` / `Sequence[int]` / `npt.ArrayLike`. Also finishes a dedup that 74327041 got wrong: the `coeffs has length ...` check appears four times in `lindblad.rs`, and the replacement there landed on the first occurrence rather than the intended one in `pc_step_orbit_rep`. All four now call `check_coeffs_len`. Messages were identical, so no behavior changed either way. 288 Python tests pass (three new: the same-object rejection and num_threads=1/2 result-equivalence), `cargo test --workspace` green, clippy/ruff/ty clean. Co-Authored-By: Claude Opus 5 (1M context) --- crates/ppvm-lindblad/src/step.rs | 16 +++++++++++ crates/ppvm-python-native/src/interface.rs | 11 +++++++- crates/ppvm-python-native/src/lindblad.rs | 28 ++++--------------- ppvm-python/src/ppvm/_core.pyi | 1 + ppvm-python/src/ppvm/lindblad.py | 11 ++++++-- ppvm-python/src/ppvm/paulisum.py | 9 ++++-- .../test/lindblad/test_pc_step_orbit_rep.py | 22 +++++++++++++++ ppvm-python/test/test_momentum_merge.py | 11 ++++++++ 8 files changed, 82 insertions(+), 27 deletions(-) diff --git a/crates/ppvm-lindblad/src/step.rs b/crates/ppvm-lindblad/src/step.rs index d0e380a0e..7a168a77e 100644 --- a/crates/ppvm-lindblad/src/step.rs +++ b/crates/ppvm-lindblad/src/step.rs @@ -198,6 +198,8 @@ impl LindbladSpec { /// representatives. If not, call /// [`canonicalize_basis_to_rep`](crate::canonicalize_basis_to_rep) /// first. + /// + /// Honours `cfg.num_threads` the same way [`Self::pc_step`] does. pub fn pc_step_orbit_rep( &self, basis: &mut Vec, @@ -206,6 +208,20 @@ impl LindbladSpec { protected: &[Word], sector: Sector<'_>, cfg: &PcStepConfig, + ) -> Result<(), Error> { + self.run_in_pool(cfg, |this| { + this.pc_step_orbit_rep_inner(basis, coeffs, dt, protected, sector, cfg) + }) + } + + fn pc_step_orbit_rep_inner( + &self, + basis: &mut Vec, + coeffs: &mut Vec>, + dt: f64, + protected: &[Word], + sector: Sector<'_>, + cfg: &PcStepConfig, ) -> Result<(), Error> { let PcStepConfig { max_basis, diff --git a/crates/ppvm-python-native/src/interface.rs b/crates/ppvm-python-native/src/interface.rs index 20450252b..6657b5451 100644 --- a/crates/ppvm-python-native/src/interface.rs +++ b/crates/ppvm-python-native/src/interface.rs @@ -102,10 +102,19 @@ macro_rules! create_interface_symmetry_methods { #[pyo3(signature = (other, group, momentum))] pub fn momentum_merge( &mut self, - mut other: pyo3::PyRefMut<'_, Self>, + other: &Bound<'_, Self>, group: &crate::symmetry::TranslationGroup, momentum: Vec, ) -> pyo3::PyResult<()> { + // `self` is already mutably borrowed, so passing the same + // object twice fails here — report that rather than letting + // PyO3's raw borrow error surface. + let mut other = other.try_borrow_mut().map_err(|_| { + pyo3::exceptions::PyValueError::new_err( + "momentum_merge: `self` and `other` must be distinct \ + PauliSum objects (got the same one twice)", + ) + })?; let n_q = group.core().n_qubits(); for (label, n) in [ ("self", self.inner.n_qubits()), diff --git a/crates/ppvm-python-native/src/lindblad.rs b/crates/ppvm-python-native/src/lindblad.rs index e56b78351..87f538cf6 100644 --- a/crates/ppvm-python-native/src/lindblad.rs +++ b/crates/ppvm-python-native/src/lindblad.rs @@ -202,13 +202,7 @@ impl LindbladSpec { let mut basis_words = decode_basis(&basis_view, n_q)?; assert_basis_unique(&basis_words)?; let mut coeffs_vec = coeffs.as_slice()?.to_vec(); - if coeffs_vec.len() != basis_words.len() { - return Err(PyValueError::new_err(format!( - "coeffs has length {} but basis has {} rows", - coeffs_vec.len(), - basis_words.len() - ))); - } + check_coeffs_len(coeffs_vec.len(), basis_words.len())?; let protected_words: Vec = if let Some(ref p) = protected { decode_basis(&p.as_array(), n_q)? } else { @@ -264,13 +258,7 @@ impl LindbladSpec { let mut basis_words = decode_basis(&basis_view, n_q)?; assert_basis_unique(&basis_words)?; let mut coeffs_vec = coeffs.as_slice()?.to_vec(); - if coeffs_vec.len() != basis_words.len() { - return Err(PyValueError::new_err(format!( - "coeffs has length {} but basis has {} rows", - coeffs_vec.len(), - basis_words.len() - ))); - } + check_coeffs_len(coeffs_vec.len(), basis_words.len())?; let protected_words: Vec = if let Some(ref p) = protected { decode_basis(&p.as_array(), n_q)? } else { @@ -334,6 +322,7 @@ impl LindbladSpec { canonicalize_first = false, admit_basis = None, tau_add = None, + num_threads = None, ))] #[allow(clippy::too_many_arguments)] fn pc_step_orbit_rep<'py>( @@ -350,6 +339,7 @@ impl LindbladSpec { canonicalize_first: bool, admit_basis: Option, tau_add: Option, + num_threads: Option, ) -> PyResult> { use num::Complex; use ppvm_lindblad::{Sector, canonicalize_basis_to_rep}; @@ -358,13 +348,7 @@ impl LindbladSpec { let basis_view = basis.as_array(); let mut basis_words = decode_basis(&basis_view, n_q)?; let coeffs_slice = coeffs.as_slice()?; - if coeffs_slice.len() != basis_words.len() { - return Err(PyValueError::new_err(format!( - "coeffs has length {} but basis has {} rows", - coeffs_slice.len(), - basis_words.len() - ))); - } + check_coeffs_len(coeffs_slice.len(), basis_words.len())?; let mut coeffs_vec: Vec> = coeffs_slice .iter() .map(|c| Complex::new(c.re, c.im)) @@ -391,7 +375,7 @@ impl LindbladSpec { admit_basis, drop_tol, tau_add, - num_threads: None, + num_threads, }, ) .map_err(map_err)?; diff --git a/ppvm-python/src/ppvm/_core.pyi b/ppvm-python/src/ppvm/_core.pyi index b0350fc55..6ea1696ca 100644 --- a/ppvm-python/src/ppvm/_core.pyi +++ b/ppvm-python/src/ppvm/_core.pyi @@ -413,6 +413,7 @@ class LindbladSpec: canonicalize_first: bool = False, admit_basis: int | None = None, tau_add: float | None = None, + num_threads: int | None = None, ) -> tuple[np.ndarray, np.ndarray]: ... def generator(self, basis: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]: ... diff --git a/ppvm-python/src/ppvm/lindblad.py b/ppvm-python/src/ppvm/lindblad.py index 82dfcd769..c2538ab5f 100644 --- a/ppvm-python/src/ppvm/lindblad.py +++ b/ppvm-python/src/ppvm/lindblad.py @@ -42,7 +42,9 @@ from typing import Union import numpy as np +import numpy.typing as npt +from . import _core from ._core import LindbladSpec as _LindbladSpec _PAULI_CODE = {"I": 0, "X": 1, "Z": 2, "Y": 3} @@ -297,13 +299,14 @@ def pc_step_orbit_rep( coeffs: np.ndarray, dt: float, max_basis: int, - group, - momentum: np.ndarray, + group: _core.TranslationGroup, + momentum: npt.ArrayLike, drop_tol: float = 1e-12, protected_arr: np.ndarray | None = None, canonicalize_first: bool = False, admit_basis: int | None = None, tau_add: float | None = None, + num_threads: int | None = None, ) -> tuple[np.ndarray, np.ndarray]: """Per-step orbit-representative pc evolution. @@ -330,6 +333,9 @@ def pc_step_orbit_rep( ``basis_arr`` is assumed to contain canonical reps only. Pass ``canonicalize_first=True`` to rewrite each row to its canonical rep on entry (coefficients unchanged). + + ``num_threads``, when set, pins this call to a freshly-built rayon + pool of that size, exactly as for `pc_step_arr`. """ n = self.n_qubits if protected_arr is None: @@ -346,6 +352,7 @@ def pc_step_orbit_rep( bool(canonicalize_first), None if admit_basis is None else int(admit_basis), None if tau_add is None else float(tau_add), + None if num_threads is None else int(num_threads), ) def pc_step( diff --git a/ppvm-python/src/ppvm/paulisum.py b/ppvm-python/src/ppvm/paulisum.py index f71c6d846..7d301fdd0 100644 --- a/ppvm-python/src/ppvm/paulisum.py +++ b/ppvm-python/src/ppvm/paulisum.py @@ -386,7 +386,7 @@ def trace(self, pattern: str) -> float: """ return self._interface.trace(pattern) - def symmetry_merge(self, group) -> None: + def symmetry_merge(self, group: _core.TranslationGroup) -> None: """Merge entries into orbit-representative form under a translation group. Each Pauli word in the sum is replaced by its canonical (lex-min) @@ -407,7 +407,12 @@ def symmetry_merge(self, group) -> None: """ self._interface.symmetry_merge(group) - def momentum_merge(self, other: "PauliSum", group, momentum) -> None: + def momentum_merge( + self, + other: "PauliSum", + group: _core.TranslationGroup, + momentum: Sequence[int], + ) -> None: """Phase-aware (momentum-sector) merge for a complex operator stored as a *real pair*: ``self`` is the real part and ``other`` the imaginary part of ``O = self + i·other``. Both are overwritten in diff --git a/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py b/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py index 36dddaf1b..b4a014d15 100644 --- a/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py +++ b/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py @@ -199,6 +199,28 @@ def test_protected_reps_are_never_dropped(): assert keep <= {string(w) for w in basis} +@pytest.mark.parametrize("num_threads", [1, 2]) +def test_num_threads_does_not_change_the_result(num_threads): + """``num_threads`` pins the call to a fresh rayon pool — same result, and + (unlike before) it is no longer silently ignored on this path.""" + n = 4 + dt = 0.03 + op = xy_chain_pbc(n, gamma=0.2) + group = TranslationGroup.chain_1d(n) + k_arr = momentum(1) + seed_basis, seed_coeffs = z_momentum_seed(n, 1) + basis, coeffs = canonicalize_basis_arr_complex(seed_basis, seed_coeffs, group, k_arr) + + ref_b, ref_c = op.pc_step_orbit_rep(basis, coeffs, dt, 10_000_000, group, k_arr) + got_b, got_c = op.pc_step_orbit_rep( + basis, coeffs, dt, 10_000_000, group, k_arr, num_threads=num_threads + ) + ref, got = to_dict(ref_b, ref_c), to_dict(got_b, got_c) + assert ref.keys() == got.keys() + for w in ref: + assert abs(ref[w] - got[w]) < 1e-12 + + def test_pc_step_orbit_rep_validates_inputs(): n = 3 op = xy_chain_pbc(n, gamma=0.0) diff --git a/ppvm-python/test/test_momentum_merge.py b/ppvm-python/test/test_momentum_merge.py index 217863299..5f8c02a70 100644 --- a/ppvm-python/test/test_momentum_merge.py +++ b/ppvm-python/test/test_momentum_merge.py @@ -115,6 +115,17 @@ def test_momentum_merge_idempotent_on_stabilized_orbit(word): assert max(abs(once.get(x, 0j) - twice.get(x, 0j)) for x in keys) < 1e-12 +def test_momentum_merge_rejects_the_same_object_twice(): + """``self`` and ``other`` hold the real and imaginary parts, so they must + be distinct; passing one object twice gets a message saying so rather + than a raw borrow error.""" + n = 4 + g = TranslationGroup.chain_1d(n) + PA, _ = _seed_pair(n, 1) + with pytest.raises(ValueError, match="must be distinct PauliSum objects"): + PA.momentum_merge(PA, g, [1]) + + def test_momentum_merge_projects_out_other_sectors(): """Merging a pure sector-k operator in sector k' != k gives ~zero.""" n = 4 From d8313e0ebf115b2f367c3884efa60e313e8138ba Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 2 Sep 2026 15:26:28 +0200 Subject: [PATCH 13/13] fix(symmetry): orbit-rep evolution on stabilized orbits; validate group inputs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two fixes from a review of the orbit-rep CTPP path. **Stabilized orbits.** The phase-aware action sums characters over the output orbit's *distinct* members, which makes it the orbit-rep generator in the summing convention `ĉ_r = |orbit_r| · c_r`. The public API carries averaged coefficients (`c_r` = the plain coefficient of the rep word, as `canonicalize_pauli_sum_complex` produces), so every matrix entry needs the similarity factor `|orbit_in| / |orbit_out|` — which is 1 only when both orbits are free. Without it, evolving `ZIZI + IZIZ` on a 4-site chain came out a factor 2 too large on every coefficient. Reps whose stabilizer is incompatible with the sector are now dropped too, matching what the reference projection does. `TranslationGroup::canonicalize_in_sector` computes the rep, the shift counter, the distinct-orbit size and sector compatibility from ONE orbit traversal (orbit-stabilizer, counting stabilizer elements during the lex-min walk), so the hot action loop pays one extra word compare. **FFI validation.** `pc_step_orbit_rep` never checked the group's qubit count against the spec's, and ran `canonicalize_first` — documented as not deduplicating — without a uniqueness check, so a release wheel silently collapsed same-orbit rows onto one CSC index. The lattice constructors and the unchecked half of `from_generators` (zero order, inexact order, non-commuting generators) aborted through `assert!`. All now raise `ValueError`. The order and commutation checks move into a fallible core constructor, `TranslationGroup::try_from_generators` -> `GroupError`, rather than being duplicated in the binding; `from_generators` panics with the same messages as before. `LossyPauliSum` gained explicit `symmetry_merge` / `momentum_merge` overrides raising `NotImplementedError`, and the `.pyi` stubs moved onto a non-loss base so type checkers reject the call. Co-Authored-By: Claude Opus 5 (1M context) --- crates/ppvm-lindblad/src/basis.rs | 22 ++- crates/ppvm-lindblad/src/mf_expm.rs | 31 ++- crates/ppvm-lindblad/src/sector.rs | 47 ++++- crates/ppvm-lindblad/src/tests.rs | 130 ++++++++----- crates/ppvm-pauli-sum/src/symmetry/group.rs | 179 ++++++++++++++---- crates/ppvm-pauli-sum/src/symmetry/mod.rs | 2 +- .../ppvm-pauli-sum/src/symmetry/momentum.rs | 48 +++++ crates/ppvm-pauli-sum/src/symmetry/tests.rs | 117 ++++++++++++ crates/ppvm-python-native/src/lindblad.rs | 9 +- crates/ppvm-python-native/src/pauli_arr.rs | 13 ++ crates/ppvm-python-native/src/symmetry.rs | 90 ++++----- ppvm-python/src/ppvm/_core.pyi | 39 ++-- ppvm-python/src/ppvm/lindblad.py | 8 + ppvm-python/src/ppvm/paulisum.py | 30 +++ .../test/lindblad/test_pc_step_orbit_rep.py | 75 ++++++++ ppvm-python/test/test_symmetry_merge.py | 42 +++- 16 files changed, 720 insertions(+), 162 deletions(-) diff --git a/crates/ppvm-lindblad/src/basis.rs b/crates/ppvm-lindblad/src/basis.rs index 941c67b6c..153780213 100644 --- a/crates/ppvm-lindblad/src/basis.rs +++ b/crates/ppvm-lindblad/src/basis.rs @@ -241,9 +241,15 @@ impl LindbladSpec { /// /// 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)` via [`Sector::canonicalize_phase`]. + /// 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`. + /// `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. /// @@ -293,12 +299,20 @@ impl LindbladSpec { |(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 (r_q, phase) = sector.canonicalize_phase(q); + let Some((r_q, phase, orbit_out)) = sector.canonicalize_phase(q) else { + continue; + }; if !in_basis.contains(&r_q) && !protected_set.contains(&r_q) { - out.push((r_q, phase * *v * c_r)); + let rate = phase * *v * c_r * (orbit_in as f64 / orbit_out as f64); + out.push((r_q, rate)); } } out diff --git a/crates/ppvm-lindblad/src/mf_expm.rs b/crates/ppvm-lindblad/src/mf_expm.rs index 491f686d5..9437396f8 100644 --- a/crates/ppvm-lindblad/src/mf_expm.rs +++ b/crates/ppvm-lindblad/src/mf_expm.rs @@ -92,11 +92,21 @@ fn build_mf_cols( /// generator `M` at momentum `sector`, plus the `(m, s)`/`μ` selection data /// — from ONE action pass over the basis. /// -/// `cols[c]` holds `(row, χ_k(g_{cnt_q}) · v_q)` for every action output -/// Pauli `q` of `L*(basis[c])` whose orbit rep `r_q` is in `basis` at index -/// `row`; outputs whose rep is out of basis are dropped. This is the -/// expensive part of the orbit-rep dynamics (`compute_action_terms`, -/// [`Sector::canonicalize_phase`]). +/// `cols[c]` holds `(row, χ_k(g_{cnt_q}) · v_q · |orbit_c| / |orbit_row|)` +/// for every action output Pauli `q` of `L*(basis[c])` whose orbit rep +/// `r_q` is in `basis` at index `row`; outputs whose rep is out of basis +/// are dropped. This is the expensive part of the orbit-rep dynamics +/// (`compute_action_terms`, [`Sector::canonicalize_phase`]). +/// +/// The character-weighted sum runs over the *output* orbit's distinct +/// members, which makes it the generator in the **summing** convention +/// `ĉ_r = |orbit_r| · c_r`. Coefficients here are in the *averaged* +/// convention (`c_r` = the plain coefficient of the rep word, what +/// `canonicalize_pauli_sum_complex` produces), so each entry carries the +/// similarity factor `|orbit_c| / |orbit_row|` that converts between +/// them. It is 1 exactly when both orbits are free — hence the factor is +/// invisible until an orbit has a non-trivial stabilizer, and cannot be +/// hoisted out as a global `|G|`. /// /// Unlike [`build_mf_cols`], `per_col[c].0` sums only the retained /// in-basis entries — the exact column 1-norm of the restricted `M`, not an @@ -124,14 +134,21 @@ fn build_orbit_rep_cols( ) }, |(s1, s2, lm), (c, r)| { + // A rep that cannot carry the sector has coefficient zero + // identically, so its column is empty. + let Some(orbit_in) = sector.orbit_size(r) else { + return (Vec::new(), (0.0, Complex::new(0.0, 0.0))); + }; let terms = spec.compute_action_terms(r, s1, s2, lm); let mut out = Vec::with_capacity(terms.len()); let mut raw = 0.0; let mut diag = Complex::new(0.0, 0.0); for (q, v) in terms.iter() { - let (r_q, phase) = sector.canonicalize_phase(q); + let Some((r_q, phase, orbit_out)) = sector.canonicalize_phase(q) else { + continue; + }; if let Some(&row) = index.get(&r_q) { - let val = phase * *v; + let val = phase * *v * (orbit_in as f64 / orbit_out as f64); raw += val.norm(); if row as usize == c { diag += val; diff --git a/crates/ppvm-lindblad/src/sector.rs b/crates/ppvm-lindblad/src/sector.rs index ab588f7aa..b76ad2687 100644 --- a/crates/ppvm-lindblad/src/sector.rs +++ b/crates/ppvm-lindblad/src/sector.rs @@ -9,9 +9,17 @@ //! representatives and the coefficients are complex (one per rep). The //! dynamics `L*` is computed with **phase-aware action** — for each //! output Pauli `q`, we canonicalize `q` to its orbit rep `r_q` with -//! shift counter `cnt_q`, and accumulate `χ_k(g_{cnt_q}) · v · c_r` -//! (where `v` is the matrix element of `L*` between input rep `r` and -//! output `q`). [`Sector::canonicalize_phase`] is that step. +//! shift counter `cnt_q`, and accumulate +//! `χ_k(g_{cnt_q}) · v · c_r · |orbit_r| / |orbit_{r_q}|` (where `v` is +//! the matrix element of `L*` between input rep `r` and output `q`). +//! [`Sector::canonicalize_phase`] is that step. +//! +//! Coefficients are in the **averaged** convention: `c_r` is the plain +//! coefficient of the rep word, as produced by +//! `canonicalize_pauli_sum_complex`. The character-weighted action is +//! naturally the generator in the *summing* convention +//! `ĉ_r = |orbit_r| · c_r`, which is where the orbit-size ratio comes +//! from; it is 1 whenever both orbits are free. //! //! The orbit-rep basis is ~`|G|`× smaller than the full-basis //! representation, throughout the entire evolution. @@ -51,13 +59,36 @@ impl<'a> Sector<'a> { /// Canonicalize `q` to its orbit representative `r_q` and return it /// alongside the character phase `χ_k(g_{cnt_q})` of the group - /// element that maps `q` to `r_q`. The phase weights the matrix - /// element of `L*` when it is accumulated onto `r_q`. + /// element that maps `q` to `r_q`, and the number of **distinct** + /// members of that orbit. The phase weights the matrix element of + /// `L*` when it is accumulated onto `r_q`; the orbit size converts + /// between the two coefficient conventions (see + /// [`Self::orbit_size`]). + /// + /// `None` when `q`'s orbit cannot carry this sector (its stabilizer + /// is incompatible with `k`): the coefficient of such a rep is + /// identically zero, so the term is dropped. #[inline] - pub fn canonicalize_phase(&self, q: &Word) -> (Word, Complex) { - let (rep, counter) = self.group.canonicalize_with_shift(q); + pub fn canonicalize_phase(&self, q: &Word) -> Option<(Word, Complex, usize)> { + let (rep, counter, orbit_size) = self.group.canonicalize_in_sector(q, self.k_modes)?; let phase = self.group.character(self.k_modes, &counter); - (rep, phase) + Some((rep, phase, orbit_size)) + } + + /// Number of **distinct** members of `w`'s translation orbit, or + /// `None` if the orbit cannot carry this sector. + /// + /// This is the factor between the two orbit-rep coefficient + /// conventions: the *averaged* one, in which `c_r` is the plain + /// coefficient of the rep word (what `canonicalize_pauli_sum_complex` + /// and this crate's public orbit-rep API use), and the *summing* one + /// `ĉ_r = |orbit_r| · c_r` (what `momentum_merge_pauli_sum_pair` + /// uses). It is `|G|` only for free orbits. + #[inline] + pub fn orbit_size(&self, w: &Word) -> Option { + self.group + .canonicalize_in_sector(w, self.k_modes) + .map(|(_, _, orbit_size)| orbit_size) } } diff --git a/crates/ppvm-lindblad/src/tests.rs b/crates/ppvm-lindblad/src/tests.rs index 5687c0b31..85526d37f 100644 --- a/crates/ppvm-lindblad/src/tests.rs +++ b/crates/ppvm-lindblad/src/tests.rs @@ -94,68 +94,61 @@ fn word_codec_roundtrip() { assert_eq!(out.as_slice(), &codes); } -/// Per-step orbit-rep evolution gives the SAME final orbit-rep -/// state as full-basis complex evolution followed by a single -/// projection at the end. Validates that the phase-aware complex -/// action machinery is consistent with the full-basis reference. -#[test] -fn pc_step_orbit_rep_matches_full_basis_projection() { - use std::f64::consts::PI; - - use ppvm_pauli_sum::symmetry::canonicalize_pauli_sum_complex; - let n = 4usize; - let dt = 0.01f64; - let n_steps = 3usize; +/// Translation-invariant XY chain with PBC on `n` sites, no dissipation. +fn xy_chain_pbc(n: usize) -> Vec<(String, f64)> { let mut h_terms: Vec<(String, f64)> = Vec::new(); for j in 0..n { let nxt = (j + 1) % n; - for op in ["X", "Y"] { + for op in ['X', 'Y'] { let mut s = vec!['I'; n]; - s[j] = op.chars().next().unwrap(); - s[nxt] = op.chars().next().unwrap(); + s[j] = op; + s[nxt] = op; h_terms.push((s.into_iter().collect(), 1.0)); } } - let spec = LindbladSpec::new(n, &h_terms, &[]).unwrap(); - let group = ppvm_pauli_sum::symmetry::TranslationGroup::chain_1d(n); - let k_mode: i32 = 1; - let k = vec![k_mode]; + h_terms +} - // Build the k=1 eigenstate in FULL basis form. - let basis_full: Vec = (0..n) - .map(|j| { - let mut s = vec!['I'; n]; - s[j] = 'Z'; - let (w, _) = parse_pauli_string(&s.into_iter().collect::(), n).unwrap(); - w - }) - .collect(); - let coeffs_full: Vec> = (0..n as i32) - .map(|a| Complex::from_polar(1.0, -2.0 * PI * (k_mode as f64) * (a as f64) / (n as f64))) +/// Per-step orbit-rep evolution must give the SAME final orbit-rep state +/// as full-basis complex evolution followed by a single projection at the +/// end (the projection theorem), for a `seed` that is a momentum-`k` +/// eigenstate in full-basis form. +/// +/// Both sides run untruncated: `pc_step_complex_full` admits every +/// leakage string, and the orbit-rep side gets a huge `max_basis`, so the +/// only remaining difference would be a bug in the phase-aware action. +fn assert_orbit_rep_matches_projection( + n: usize, + h_terms: &[(String, f64)], + seed: &[(&str, Complex)], + k: &[i32], + dt: f64, + n_steps: usize, +) { + use ppvm_pauli_sum::symmetry::canonicalize_pauli_sum_complex; + + let spec = LindbladSpec::new(n, h_terms, &[]).unwrap(); + let group = ppvm_pauli_sum::symmetry::TranslationGroup::chain_1d(n); + let basis_full: Vec = seed + .iter() + .map(|(s, _)| parse_pauli_string(s, n).unwrap().0) .collect(); + let coeffs_full: Vec> = seed.iter().map(|(_, c)| *c).collect(); - // ----- Full-basis path ----- + // ----- Full-basis path, projected once at the end ----- let mut bf = basis_full.clone(); let mut cf = coeffs_full.clone(); let protected: Vec = Vec::new(); for _ in 0..n_steps { - // Full enrichment (tau_add = 0.0 adds every leakage string): - // for a momentum eigenstate the leakage is pure-sector, so the - // full-basis and orbit-rep paths build corresponding bases and - // the projection theorem gives an exact match. The orbit-rep - // side uses a large max_basis so its rank cap never binds. pc_step_complex_full(&spec, &mut bf, &mut cf, dt); } - // Project at the end. - canonicalize_pauli_sum_complex(&mut bf, &mut cf, &group, &k); + canonicalize_pauli_sum_complex(&mut bf, &mut cf, &group, k); - // ----- Orbit-rep path ----- - // Initial orbit-rep form: project the full-basis input. + // ----- Orbit-rep path: project the seed, then evolve in rep form ----- let mut br = basis_full.clone(); let mut cr = coeffs_full.clone(); - canonicalize_pauli_sum_complex(&mut br, &mut cr, &group, &k); - // Evolve in orbit-rep form (max_basis large ⇒ full enrichment). - let sector = Sector::new(&group, &k); + canonicalize_pauli_sum_complex(&mut br, &mut cr, &group, k); + let sector = Sector::new(&group, k); for _ in 0..n_steps { spec.pc_step_orbit_rep( &mut br, @@ -171,7 +164,6 @@ fn pc_step_orbit_rep_matches_full_basis_projection() { .unwrap(); } - // Compare. let mf: FxHashMap> = bf.into_iter().zip(cf).collect(); let mr: FxHashMap> = br.into_iter().zip(cr).collect(); assert_eq!( @@ -186,7 +178,7 @@ fn pc_step_orbit_rep_matches_full_basis_projection() { let cf_val = mf .get(w) .copied() - .unwrap_or_else(|| panic!("rep {:?} in orbit-rep but not in full-basis", w)); + .unwrap_or_else(|| panic!("rep {w} in orbit-rep but not in full-basis")); max_diff = max_diff.max((cm - cf_val).norm()); } assert!( @@ -195,6 +187,54 @@ fn pc_step_orbit_rep_matches_full_basis_projection() { ); } +/// Validates the phase-aware complex action against the full-basis +/// reference on a `k=1` seed whose orbits are all free. +#[test] +fn pc_step_orbit_rep_matches_full_basis_projection() { + use std::f64::consts::PI; + + let n = 4usize; + let k_mode = 1i32; + // `O_k = Σ_a e^{-2πi k a / n} Z_a`, a k=1 momentum eigenstate. + let words: Vec = (0..n) + .map(|j| { + let mut s = vec!['I'; n]; + s[j] = 'Z'; + s.into_iter().collect() + }) + .collect(); + let seed: Vec<(&str, Complex)> = words + .iter() + .enumerate() + .map(|(a, s)| { + let phase = -2.0 * PI * (k_mode as f64) * (a as f64) / (n as f64); + (s.as_str(), Complex::from_polar(1.0, phase)) + }) + .collect(); + + assert_orbit_rep_matches_projection(n, &xy_chain_pbc(n), &seed, &[k_mode], 0.01, 3); +} + +/// Same projection-theorem check on a seed living on a **stabilized** +/// orbit: `ZIZI + IZIZ` has period 2 on a 4-site chain, so its orbit has +/// 2 distinct members, not 4. +/// +/// Regression test: the phase-aware action is the orbit-rep generator in +/// the *summing* convention, and converting it to the *averaged* +/// convention that `canonicalize_pauli_sum_complex` uses takes a per- +/// orbit-pair `|orbit_in| / |orbit_out|` factor — which is 1 only when +/// both orbits are free. Without that factor this evolves `ZIZI + IZIZ` +/// with every coefficient off by exactly 2. +#[test] +fn pc_step_orbit_rep_handles_stabilized_orbits() { + let n = 4usize; + let seed = [ + ("ZIZI", Complex::new(1.0, 0.0)), + ("IZIZ", Complex::new(1.0, 0.0)), + ]; + assert_orbit_rep_matches_projection(n, &xy_chain_pbc(n), &seed, &[0], 0.01, 3); +} + /// The full-space complex step at momentum k=0 must reproduce the real /// pc_step on the same trajectory exactly. #[test] diff --git a/crates/ppvm-pauli-sum/src/symmetry/group.rs b/crates/ppvm-pauli-sum/src/symmetry/group.rs index cda4ae7c7..73deab159 100644 --- a/crates/ppvm-pauli-sum/src/symmetry/group.rs +++ b/crates/ppvm-pauli-sum/src/symmetry/group.rs @@ -62,6 +62,92 @@ pub(super) fn validate_site_count(n: usize, context: &str) { .unwrap_or_else(|_| panic!("{context}: site count {n} exceeds the u32-addressable range")); } +/// A precondition violation in [`TranslationGroup::try_from_generators`]. +/// +/// Every variant is caller-supplied-input error, and its [`Display`] +/// text is exactly what [`TranslationGroup::from_generators`] panics +/// with. Arithmetic overflow in the group order or character phase +/// modulus is NOT covered — that needs generator orders in the billions +/// and still panics. +/// +/// [`Display`]: std::fmt::Display +#[derive(Debug, Clone, PartialEq, Eq)] +pub enum GroupError { + /// `perms` and `orders` describe different numbers of generators. + LengthMismatch { perms: usize, orders: usize }, + /// A generator's permutation is not `n_qubits` long. + PermutationLength { + generator: usize, + len: usize, + n_qubits: usize, + }, + /// A generator maps a qubit outside `0..n_qubits`. + TargetOutOfRange { + generator: usize, + target: u32, + n_qubits: usize, + }, + /// A generator maps two qubits to the same position. + DuplicateTarget { generator: usize, target: u32 }, + /// A generator declares cyclic order zero. + ZeroOrder { generator: usize }, + /// A generator's declared order is not its exact cyclic order. + OrderMismatch { + generator: usize, + declared: u32, + exact: u32, + }, + /// Two generators do not commute, so they generate no abelian group. + NonCommuting { left: usize, right: usize }, +} + +impl std::fmt::Display for GroupError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::LengthMismatch { perms, orders } => write!( + f, + "perms ({perms} generators) and orders ({orders}) must have the same length" + ), + Self::PermutationLength { + generator, + len, + n_qubits, + } => write!( + f, + "generator {generator}: permutation length {len} != n_qubits {n_qubits}" + ), + Self::TargetOutOfRange { + generator, + target, + n_qubits, + } => write!( + f, + "generator {generator}: target {target} out of range [0, {n_qubits})" + ), + Self::DuplicateTarget { generator, target } => write!( + f, + "generator {generator}: not a permutation (duplicate target {target})" + ), + Self::ZeroOrder { generator } => { + write!(f, "generator {generator} order must be nonzero") + } + Self::OrderMismatch { + generator, + declared, + exact, + } => write!( + f, + "generator {generator} declared order {declared} != exact permutation order {exact}" + ), + Self::NonCommuting { left, right } => { + write!(f, "generators {left} and {right} do not commute") + } + } + } +} + +impl std::error::Error for GroupError {} + /// A finite abelian symmetry group acting on qubit positions by /// permutations. /// @@ -94,61 +180,90 @@ pub struct TranslationGroup { } impl TranslationGroup { - /// Construct from explicit generator permutations and orders. + /// Construct from explicit generator permutations and orders, + /// panicking on any precondition violation. /// /// Each `perm` must be a permutation of `0..n_qubits`. Each `order` /// must be the permutation's exact cyclic order, not merely a /// multiple for which `perm^order == identity`. Generators must /// commute, but their combined action may still have a kernel. + /// + /// Use [`Self::try_from_generators`] when the generators come from + /// outside the program (an FFI boundary, a config file) and a + /// precondition violation should be reported rather than abort. pub fn from_generators(n_qubits: usize, perms: Vec>, orders: Vec) -> Self { - assert_eq!(perms.len(), orders.len(), "perms and orders must match"); - for (g, perm) in perms.iter().enumerate() { - assert_eq!( - perm.len(), - n_qubits, - "generator {g} permutation has length {} != n_qubits {n_qubits}", - perm.len() - ); + Self::try_from_generators(n_qubits, perms, orders) + .unwrap_or_else(|err| panic!("TranslationGroup::from_generators: {err}")) + } + + /// Fallible [`Self::from_generators`]: validates every precondition + /// on the caller-supplied generators and reports the first violation + /// as a [`GroupError`] instead of panicking. + pub fn try_from_generators( + n_qubits: usize, + perms: Vec>, + orders: Vec, + ) -> Result { + if perms.len() != orders.len() { + return Err(GroupError::LengthMismatch { + perms: perms.len(), + orders: orders.len(), + }); + } + for (generator, perm) in perms.iter().enumerate() { + if perm.len() != n_qubits { + return Err(GroupError::PermutationLength { + generator, + len: perm.len(), + n_qubits, + }); + } let mut seen = vec![false; n_qubits]; - for &p in perm { - assert!( - (p as usize) < n_qubits, - "generator {g} maps to out-of-range position {p}" - ); - assert!( - !seen[p as usize], - "generator {g} is not a permutation (duplicate target {p})" - ); - seen[p as usize] = true; + for &target in perm { + if target as usize >= n_qubits { + return Err(GroupError::TargetOutOfRange { + generator, + target, + n_qubits, + }); + } + if seen[target as usize] { + return Err(GroupError::DuplicateTarget { generator, target }); + } + seen[target as usize] = true; } } - for (g, &declared) in orders.iter().enumerate() { - assert!(declared != 0, "generator {g} order must be nonzero"); - let exact = permutation_order(&perms[g], g); - assert_eq!( - declared, exact, - "generator {g} declared order {declared} != exact permutation order {exact}", - ); + for (generator, &declared) in orders.iter().enumerate() { + if declared == 0 { + return Err(GroupError::ZeroOrder { generator }); + } + let exact = permutation_order(&perms[generator], generator); + if declared != exact { + return Err(GroupError::OrderMismatch { + generator, + declared, + exact, + }); + } } for left in 0..perms.len() { for right in left + 1..perms.len() { - assert!( - permutations_commute(&perms[left], &perms[right]), - "generators {left} and {right} do not commute", - ); + if !permutations_commute(&perms[left], &perms[right]) { + return Err(GroupError::NonCommuting { left, right }); + } } } let order = checked_group_order(&orders); let phase_modulus = orders.iter().fold(1usize, |acc, &value| { checked_lcm(acc, value as usize, "character phase modulus") }); - Self { + Ok(Self { n_qubits, perms, orders, order, phase_modulus, - } + }) } /// 1D chain of `n` sites with periodic boundary conditions. diff --git a/crates/ppvm-pauli-sum/src/symmetry/mod.rs b/crates/ppvm-pauli-sum/src/symmetry/mod.rs index 22fdbb291..b8ee315e1 100644 --- a/crates/ppvm-pauli-sum/src/symmetry/mod.rs +++ b/crates/ppvm-pauli-sum/src/symmetry/mod.rs @@ -72,7 +72,7 @@ mod group; mod merge; mod momentum; -pub use group::TranslationGroup; +pub use group::{GroupError, TranslationGroup}; pub use merge::{canonicalize_pauli_sum, symmetry_merge_pauli_sum}; pub use momentum::{ SectorCheckError, canonicalize_pauli_sum_complex, check_momentum_sector, diff --git a/crates/ppvm-pauli-sum/src/symmetry/momentum.rs b/crates/ppvm-pauli-sum/src/symmetry/momentum.rs index 19563f27e..501cbceaa 100644 --- a/crates/ppvm-pauli-sum/src/symmetry/momentum.rs +++ b/crates/ppvm-pauli-sum/src/symmetry/momentum.rs @@ -51,6 +51,54 @@ impl TranslationGroup { let phase = 2.0 * PI * numerator as f64 / self.phase_modulus() as f64; Complex::from_polar(1.0, phase) } + + /// Everything the phase-aware routines need about `w`'s orbit in + /// momentum sector `k_modes`, from ONE orbit traversal: the lex-min + /// representative `r`, the mixed-radix counter of the group element + /// mapping `r` to `w` (as [`Self::canonicalize_with_shift`]), and the + /// number of **distinct** orbit members `|orbit|`. + /// + /// Returns `None` when the orbit's stabilizer is incompatible with + /// `k_modes` — i.e. some `s` with `s·w = w` has `χ_k(s) ≠ 1`. Such an + /// orbit cannot carry this sector: its momentum projection is + /// identically zero, and the rep coefficient a single traversal would + /// report depends on which counter the traversal happens to pick. + /// + /// `|orbit| = |G| / |stabilizer|` (orbit-stabilizer), and equals + /// `|G|` only for free orbits. + /// + /// Same `O(|G| × n_qubits)` cost as [`Self::canonicalize_with_shift`]. + pub fn canonicalize_in_sector( + &self, + w: &PauliWord, + k_modes: &[i32], + ) -> Option<(PauliWord, Vec, usize)> + where + A: PauliStorage, + S: BuildHasher + Clone + Default + HashFinalize, + { + let mut best: Option<(PauliWord, Vec)> = None; + let mut stabilizer = 0usize; + for (candidate, counter) in self.orbit_with_counters(w) { + if candidate == *w { + if self.character_numerator(k_modes, &counter) != 0 { + return None; + } + stabilizer += 1; + } + if best.as_ref().is_none_or(|(b, _)| candidate < *b) { + best = Some((candidate, counter)); + } + } + let (rep, counter_from_word) = best.expect("a finite group contains the identity element"); + let shift = (0..self.n_generators()) + .map(|g| { + let order = self.generator_order(g); + (order - counter_from_word[g]) % order + }) + .collect(); + Some((rep, shift, self.order() / stabilizer)) + } } /// Replace `(basis, complex_coeffs)` in-place with the orbit-rep form diff --git a/crates/ppvm-pauli-sum/src/symmetry/tests.rs b/crates/ppvm-pauli-sum/src/symmetry/tests.rs index 8af9a01bf..66c261d3f 100644 --- a/crates/ppvm-pauli-sum/src/symmetry/tests.rs +++ b/crates/ppvm-pauli-sum/src/symmetry/tests.rs @@ -136,6 +136,46 @@ fn canonicalize_with_shift_round_trip() { } } +#[test] +fn canonicalize_in_sector_agrees_with_canonicalize_with_shift() { + let g = TranslationGroup::chain_1d(4); + for src in ["IIXY", "IXYI", "XYII", "YIIX", "XIXI", "IIII"] { + let w = word(src); + let (rep, shift, orbit_size) = g.canonicalize_in_sector(&w, &[0]).unwrap(); + let (ref_rep, ref_shift) = g.canonicalize_with_shift(&w); + assert_eq!(rep, ref_rep, "{src}: rep"); + assert_eq!(shift, ref_shift, "{src}: shift"); + let distinct: std::collections::HashSet = g.orbit(&w).collect(); + assert_eq!(orbit_size, distinct.len(), "{src}: orbit size"); + } +} + +#[test] +fn canonicalize_in_sector_rejects_incompatible_stabilizer() { + // "XIXI" has period 2 on a 4-site chain: 2 distinct orbit members, + // stabilizer generated by T². χ_k(T²) = e^{iπk}, so the orbit + // carries the k=0 and k=2 sectors but not k=1 or k=3. + let g = TranslationGroup::chain_1d(4); + let w = word("XIXI"); + for k in [0, 2] { + let (_, _, orbit_size) = g + .canonicalize_in_sector(&w, &[k]) + .unwrap_or_else(|| panic!("k={k} must be compatible with a period-2 orbit")); + assert_eq!(orbit_size, 2, "k={k}"); + } + for k in [1, 3] { + assert!( + g.canonicalize_in_sector(&w, &[k]).is_none(), + "k={k} must be rejected on a period-2 orbit" + ); + } + // A free orbit carries every sector, with the full |G| members. + for k in 0..4 { + let (_, _, orbit_size) = g.canonicalize_in_sector(&word("XIII"), &[k]).unwrap(); + assert_eq!(orbit_size, 4, "k={k}"); + } +} + #[test] fn character_trivial_sector_is_one() { let g = TranslationGroup::chain_1d(4); @@ -626,6 +666,83 @@ fn rejects_noncommuting_generators() { TranslationGroup::from_generators(3, vec![swap_01, swap_12], vec![2, 2]); } +#[test] +fn try_from_generators_reports_every_precondition() { + use super::GroupError; + /// `(n_qubits, perms, orders, expected error)` + type Case = (usize, Vec>, Vec, GroupError); + let cases: Vec = vec![ + ( + 2, + vec![vec![1, 0]], + vec![2, 2], + GroupError::LengthMismatch { + perms: 1, + orders: 2, + }, + ), + ( + 3, + vec![vec![1, 0]], + vec![2], + GroupError::PermutationLength { + generator: 0, + len: 2, + n_qubits: 3, + }, + ), + ( + 2, + vec![vec![1, 5]], + vec![2], + GroupError::TargetOutOfRange { + generator: 0, + target: 5, + n_qubits: 2, + }, + ), + ( + 2, + vec![vec![1, 1]], + vec![2], + GroupError::DuplicateTarget { + generator: 0, + target: 1, + }, + ), + ( + 2, + vec![vec![1, 0]], + vec![0], + GroupError::ZeroOrder { generator: 0 }, + ), + ( + 2, + vec![vec![1, 0]], + vec![4], + GroupError::OrderMismatch { + generator: 0, + declared: 4, + exact: 2, + }, + ), + ( + 3, + vec![vec![1, 0, 2], vec![0, 2, 1]], + vec![2, 2], + GroupError::NonCommuting { left: 0, right: 1 }, + ), + ]; + for (n_qubits, perms, orders, expected) in cases { + let err = TranslationGroup::try_from_generators(n_qubits, perms, orders) + .expect_err("must be rejected"); + assert_eq!(err, expected); + } + // Valid input still constructs, and matches the panicking constructor. + let group = TranslationGroup::try_from_generators(4, vec![vec![1, 2, 3, 0]], vec![4]).unwrap(); + assert_eq!(group.order(), TranslationGroup::chain_1d(4).order()); +} + #[test] fn rejects_zero_lattice_dimensions() { assert!(std::panic::catch_unwind(|| TranslationGroup::chain_1d(0)).is_err()); diff --git a/crates/ppvm-python-native/src/lindblad.rs b/crates/ppvm-python-native/src/lindblad.rs index 87f538cf6..477e1e834 100644 --- a/crates/ppvm-python-native/src/lindblad.rs +++ b/crates/ppvm-python-native/src/lindblad.rs @@ -44,7 +44,9 @@ fn assert_basis_unique(basis: &[Word]) -> PyResult<()> { Ok(()) } -use crate::pauli_arr::{check_coeffs_len, check_momentum_len, decode_basis, encode_basis}; +use crate::pauli_arr::{ + check_coeffs_len, check_group_qubits, check_momentum_len, decode_basis, encode_basis, +}; /// Pack `Vec<(Word, f64)>` into the standard PyO3 return shape. fn pack_pauli_map<'py>( @@ -360,9 +362,14 @@ impl LindbladSpec { }; let k_slice = momentum.as_slice()?; check_momentum_len(k_slice.len(), group.core().n_generators())?; + check_group_qubits(n_q, group.core().n_qubits())?; if canonicalize_first { canonicalize_basis_to_rep(&mut basis_words, group.core()); } + // Canonicalization can collapse several input rows onto one rep, + // and the step indexes the basis by Pauli word — so uniqueness is + // checked after the rewrite, not before. + assert_basis_unique(&basis_words)?; self.inner .pc_step_orbit_rep( &mut basis_words, diff --git a/crates/ppvm-python-native/src/pauli_arr.rs b/crates/ppvm-python-native/src/pauli_arr.rs index a90402cc4..45580d638 100644 --- a/crates/ppvm-python-native/src/pauli_arr.rs +++ b/crates/ppvm-python-native/src/pauli_arr.rs @@ -73,6 +73,19 @@ pub(crate) fn check_group_width( Ok(()) } +/// Check that a [`crate::symmetry::TranslationGroup`] acts on the same +/// qubit count as the object being evolved. The core group routines +/// assert this on the first Pauli word they see, so without this the +/// mismatch surfaces as a panic from deep inside the step. +pub(crate) fn check_group_qubits(n_qubits: usize, group_n_qubits: usize) -> PyResult<()> { + if n_qubits != group_n_qubits { + return Err(PyValueError::new_err(format!( + "spec has {n_qubits} qubits but the TranslationGroup acts on {group_n_qubits}" + ))); + } + Ok(()) +} + /// Check that a coefficient vector has one entry per basis row. pub(crate) fn check_coeffs_len(n_coeffs: usize, n_rows: usize) -> PyResult<()> { if n_coeffs != n_rows { diff --git a/crates/ppvm-python-native/src/symmetry.rs b/crates/ppvm-python-native/src/symmetry.rs index 6c060ec51..cd3dfa32b 100644 --- a/crates/ppvm-python-native/src/symmetry.rs +++ b/crates/ppvm-python-native/src/symmetry.rs @@ -53,75 +53,75 @@ impl TranslationGroup { } } +/// Validate lattice extents before handing them to a core constructor, +/// which asserts these preconditions rather than reporting them. Each +/// extent must be positive and `u32`-addressable, and the qubit count +/// (their product) must not overflow. +fn check_lattice(dims: &[(&str, usize)]) -> PyResult<()> { + let mut n_qubits = 1usize; + for &(name, dim) in dims { + if dim == 0 { + return Err(PyValueError::new_err(format!("{name} must be positive"))); + } + n_qubits = n_qubits + .checked_mul(dim) + .ok_or_else(|| PyValueError::new_err(format!("qubit count overflows: {name}={dim}")))?; + } + u32::try_from(n_qubits - 1).map_err(|_| { + PyValueError::new_err(format!( + "qubit count {n_qubits} exceeds the u32-addressable range" + )) + })?; + Ok(()) +} + #[pymethods] impl TranslationGroup { #[staticmethod] - pub fn chain_1d(n: usize) -> Self { - Self { + pub fn chain_1d(n: usize) -> PyResult { + check_lattice(&[("n", n)])?; + Ok(Self { inner: core_sym::TranslationGroup::chain_1d(n), - } + }) } #[staticmethod] - pub fn torus_2d(lx: usize, ly: usize) -> Self { - Self { + pub fn torus_2d(lx: usize, ly: usize) -> PyResult { + check_lattice(&[("lx", lx), ("ly", ly)])?; + Ok(Self { inner: core_sym::TranslationGroup::torus_2d(lx, ly), - } + }) } #[staticmethod] - pub fn torus_3d(lx: usize, ly: usize, lz: usize) -> Self { - Self { + pub fn torus_3d(lx: usize, ly: usize, lz: usize) -> PyResult { + check_lattice(&[("lx", lx), ("ly", ly), ("lz", lz)])?; + Ok(Self { inner: core_sym::TranslationGroup::torus_3d(lx, ly, lz), - } + }) } #[staticmethod] - pub fn ladder(l: usize, n_legs: usize) -> Self { - Self { + pub fn ladder(l: usize, n_legs: usize) -> PyResult { + check_lattice(&[("l", l), ("n_legs", n_legs)])?; + Ok(Self { inner: core_sym::TranslationGroup::ladder(l, n_legs), - } + }) } + /// Every precondition — permutation shape and validity, exact cyclic + /// orders, pairwise commutation — is checked by the core's fallible + /// constructor, so bad generators raise `ValueError` here instead of + /// aborting. #[staticmethod] pub fn from_generators( n_qubits: usize, perms: Vec>, orders: Vec, ) -> PyResult { - if perms.len() != orders.len() { - return Err(PyValueError::new_err(format!( - "perms ({} generators) and orders ({}) must have the same length", - perms.len(), - orders.len() - ))); - } - for (g, perm) in perms.iter().enumerate() { - if perm.len() != n_qubits { - return Err(PyValueError::new_err(format!( - "generator {g}: permutation length {} != n_qubits {n_qubits}", - perm.len() - ))); - } - let mut seen = vec![false; n_qubits]; - for &p in perm { - let p = p as usize; - if p >= n_qubits { - return Err(PyValueError::new_err(format!( - "generator {g}: target {p} out of range [0, {n_qubits})" - ))); - } - if seen[p] { - return Err(PyValueError::new_err(format!( - "generator {g}: not a permutation (duplicate target {p})" - ))); - } - seen[p] = true; - } - } - Ok(Self { - inner: core_sym::TranslationGroup::from_generators(n_qubits, perms, orders), - }) + let inner = core_sym::TranslationGroup::try_from_generators(n_qubits, perms, orders) + .map_err(|e| PyValueError::new_err(e.to_string()))?; + Ok(Self { inner }) } /// Number of qubits this group acts on. diff --git a/ppvm-python/src/ppvm/_core.pyi b/ppvm-python/src/ppvm/_core.pyi index 6ea1696ca..e9877b7f8 100644 --- a/ppvm-python/src/ppvm/_core.pyi +++ b/ppvm-python/src/ppvm/_core.pyi @@ -62,11 +62,14 @@ class _PauliSumBase: def terms(self) -> list[tuple[str, float]]: ... def weights(self) -> list[tuple[str, int]]: ... def current_max_weight(self) -> int: ... - # Only on non-loss variants (see create_interface_symmetry_methods). + +class _PauliSumNoLossBase(_PauliSumBase): + """Methods the interface macros expand only for non-loss variants.""" + def symmetry_merge(self, group: TranslationGroup) -> None: ... def momentum_merge( self, - other: _PauliSumBase, + other: _PauliSumNoLossBase, group: TranslationGroup, momentum: list[int], ) -> None: ... @@ -78,22 +81,22 @@ class _PauliSumLossBase(_PauliSumBase): ) -> None: ... def reset_loss_channel(self, addr0: int, truncate: bool = True) -> None: ... -class PauliSumIndexMapFxHash0(_PauliSumBase): ... -class PauliSumIndexMapFxHash1(_PauliSumBase): ... -class PauliSumIndexMapFxHash2(_PauliSumBase): ... -class PauliSumIndexMapFxHash3(_PauliSumBase): ... -class PauliSumIndexMapFxHash4(_PauliSumBase): ... -class PauliSumIndexMapFxHash5(_PauliSumBase): ... -class PauliSumIndexMapFxHash6(_PauliSumBase): ... -class PauliSumIndexMapFxHash7(_PauliSumBase): ... -class PauliSumIndexMapFxHash8(_PauliSumBase): ... -class PauliSumIndexMapFxHash9(_PauliSumBase): ... -class PauliSumIndexMapFxHash10(_PauliSumBase): ... -class PauliSumIndexMapFxHash11(_PauliSumBase): ... -class PauliSumIndexMapFxHash12(_PauliSumBase): ... -class PauliSumIndexMapFxHash13(_PauliSumBase): ... -class PauliSumIndexMapFxHash14(_PauliSumBase): ... -class PauliSumIndexMapFxHash15(_PauliSumBase): ... +class PauliSumIndexMapFxHash0(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash1(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash2(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash3(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash4(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash5(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash6(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash7(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash8(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash9(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash10(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash11(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash12(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash13(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash14(_PauliSumNoLossBase): ... +class PauliSumIndexMapFxHash15(_PauliSumNoLossBase): ... class PauliSumLossIndexMapFxHash0(_PauliSumLossBase): ... class PauliSumLossIndexMapFxHash1(_PauliSumLossBase): ... class PauliSumLossIndexMapFxHash2(_PauliSumLossBase): ... diff --git a/ppvm-python/src/ppvm/lindblad.py b/ppvm-python/src/ppvm/lindblad.py index c2538ab5f..a71dce745 100644 --- a/ppvm-python/src/ppvm/lindblad.py +++ b/ppvm-python/src/ppvm/lindblad.py @@ -316,6 +316,14 @@ def pc_step_orbit_rep( is ~``|group|×`` smaller than the equivalent full-basis complex evolution, and the reduction persists across every step. + Coefficients use the same convention as + `ppvm.canonicalize_basis_arr_complex`: ``coeffs[i]`` is the plain + coefficient of the representative Pauli word itself (the + orbit-*averaged* convention, not the summing one + `PauliSum.momentum_merge` uses). Reps whose orbit cannot carry + ``momentum`` — its stabilizer has a non-trivial character — are + dropped, matching that projection. + Truncation. ``max_basis`` is a hard rank cap on the live orbit-rep basis: enrichment adds at most ``max_basis - len(basis)`` of the largest leakage reps, and the post-step basis is trimmed to the diff --git a/ppvm-python/src/ppvm/paulisum.py b/ppvm-python/src/ppvm/paulisum.py index 7d301fdd0..6d1922915 100644 --- a/ppvm-python/src/ppvm/paulisum.py +++ b/ppvm-python/src/ppvm/paulisum.py @@ -551,3 +551,33 @@ def reset_loss_channel(self, addr0: int, *, truncate: bool = True) -> None: strategy after the channel; if ``False``, defer it. """ self._interface.reset_loss_channel(addr0, truncate=truncate) + + def symmetry_merge(self, group: _core.TranslationGroup) -> None: + """Not available on `LossyPauliSum`. + + Raises: + NotImplementedError: always. Canonicalizing a lossy Pauli word + would have to permute the loss bitmap along with the Pauli + alphabet, which the Rust core does not implement. + """ + raise NotImplementedError( + "symmetry_merge is not implemented for LossyPauliSum: canonicalizing a " + "lossy Pauli word would have to permute the loss bitmap too" + ) + + def momentum_merge( + self, + other: "PauliSum", + group: _core.TranslationGroup, + momentum: Sequence[int], + ) -> None: + """Not available on `LossyPauliSum`. + + Raises: + NotImplementedError: always, for the same reason as + `symmetry_merge`. + """ + raise NotImplementedError( + "momentum_merge is not implemented for LossyPauliSum: canonicalizing a " + "lossy Pauli word would have to permute the loss bitmap too" + ) diff --git a/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py b/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py index b4a014d15..6bb8a81d2 100644 --- a/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py +++ b/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py @@ -144,6 +144,59 @@ def test_orbit_rep_matches_dense_full_space_then_project(k): assert any(abs(c) > 1e-6 for c in got.values()), "orbit-rep state decayed away" +def test_orbit_rep_handles_stabilized_orbits(): + """Same dense cross-check, seeded on a **stabilized** orbit. + + ``ZIZI + IZIZ`` has period 2 on a 4-site chain, so its orbit has 2 + distinct members rather than 4. The phase-aware action is naturally the + orbit-rep generator in the *summing* convention; converting it to the + *averaged* convention that ``canonicalize_basis_arr_complex`` returns + costs a per-orbit-pair ``|orbit_in| / |orbit_out|`` factor, which is 1 + only when both orbits are free. Without it every coefficient here comes + out exactly 2x too large. + """ + n = 4 + dt = 0.05 + n_steps = 2 + op = xy_chain_pbc(n, gamma=0.3) + group = TranslationGroup.chain_1d(n) + k_arr = momentum(0) + seed_basis, seed_coeffs = basis_arr(["ZIZI", "IZIZ"], n), np.array([1.0 + 0j, 1.0 + 0j]) + + # --- dense full-space reference --- + full = all_strings(n) + generator = np.zeros((len(full), len(full)), dtype=float) + rows, cols, vals = op.generator(full) + generator[rows, cols] = vals + index = {s: i for i, s in enumerate(full)} + v = np.zeros(len(full), dtype=complex) + for w, c in zip(seed_basis, seed_coeffs, strict=True): + v[index[string(w)]] = c + step = _dense_expm(dt * generator) + v_re, v_im = v.real.copy(), v.imag.copy() + for _ in range(n_steps): + v_re = step @ v_re + v_im = step @ v_im + expected = to_dict( + *canonicalize_basis_arr_complex(basis_arr(full, n), v_re + 1j * v_im, group, k_arr) + ) + + # --- orbit-rep evolution --- + rep_basis, rep_coeffs = canonicalize_basis_arr_complex(seed_basis, seed_coeffs, group, k_arr) + assert len(rep_coeffs) == 1, "the seed is a single orbit" + for _ in range(n_steps): + rep_basis, rep_coeffs = op.pc_step_orbit_rep( + rep_basis, rep_coeffs, dt, 10_000_000, group, k_arr, drop_tol=0.0 + ) + got = to_dict(rep_basis, rep_coeffs) + + for word in set(expected) | set(got): + e = expected.get(word, 0.0) + g = got.get(word, 0.0) + assert abs(e - g) < 1e-9, f"rep {word} dense {e} vs orbit-rep {g}" + assert any(abs(c) > 1e-6 for c in got.values()), "orbit-rep state decayed away" + + def test_canonicalize_first_accepts_non_canonical_input(): """The same physical state seeded on a non-canonical orbit member gives the same evolution once ``canonicalize_first=True`` normalizes it.""" @@ -230,6 +283,28 @@ def test_pc_step_orbit_rep_validates_inputs(): op.pc_step_orbit_rep(basis, coeffs, 0.01, 100, group, momentum(0, 0)) with pytest.raises(ValueError, match="coeffs has length 2 but basis has 3 rows"): op.pc_step_orbit_rep(basis, coeffs[:2], 0.01, 100, group, momentum(0)) + with pytest.raises(ValueError, match="spec has 3 qubits but the TranslationGroup acts on 4"): + op.pc_step_orbit_rep(basis, coeffs, 0.01, 100, TranslationGroup.chain_1d(4), momentum(0)) + + +def test_pc_step_orbit_rep_rejects_duplicate_reps(): + """The step indexes the basis by Pauli word, so duplicate rows would + silently collapse. They are rejected — including duplicates created by + ``canonicalize_first``, which does not deduplicate.""" + n = 4 + op = xy_chain_pbc(n, gamma=0.0) + group = TranslationGroup.chain_1d(n) + coeffs = np.array([1.0 + 0j, 1.0 + 0j]) + duplicate = basis_arr(["ZIZI", "ZIZI"], n) + with pytest.raises(ValueError, match="duplicate Pauli word at row 0 and row 1"): + op.pc_step_orbit_rep(duplicate, coeffs, 0.01, 100, group, momentum(0)) + # "ZIZI" and "IZIZ" are distinct words on one orbit: legal as input, + # but canonicalize_first collapses them onto the same rep. + same_orbit = basis_arr(["ZIZI", "IZIZ"], n) + with pytest.raises(ValueError, match="duplicate Pauli word at row 0 and row 1"): + op.pc_step_orbit_rep( + same_orbit, coeffs, 0.01, 100, group, momentum(0), canonicalize_first=True + ) def test_returns_complex_arrays_of_matching_shape(): diff --git a/ppvm-python/test/test_symmetry_merge.py b/ppvm-python/test/test_symmetry_merge.py index 4de560263..f91536534 100644 --- a/ppvm-python/test/test_symmetry_merge.py +++ b/ppvm-python/test/test_symmetry_merge.py @@ -12,7 +12,7 @@ import numpy as np import pytest -from ppvm import PauliSum, TranslationGroup +from ppvm import LossyPauliSum, PauliSum, TranslationGroup _CODE = {"I": 0, "X": 1, "Z": 2, "Y": 3} _CHAR = {v: k for k, v in _CODE.items()} @@ -63,13 +63,40 @@ def test_from_generators_matches_chain_1d(): ([[1, 0, 2]], [2], "permutation length"), ([[1, 0, 2, 9]], [2], "out of range"), ([[1, 1, 2, 3]], [2], "duplicate target"), + # The declared order must be the permutation's *exact* cyclic + # order, not a multiple of it and not zero. + ([[1, 2, 3, 0]], [2], "declared order 2 != exact permutation order 4"), + ([[1, 2, 3, 0]], [8], "declared order 8 != exact permutation order 4"), + ([[1, 2, 3, 0]], [0], "order must be nonzero"), + # Generators must commute: (0 1) and (1 2) do not. + ([[1, 0, 2, 3], [0, 2, 1, 3]], [2, 2], "generators 0 and 1 do not commute"), ], ) def test_from_generators_validates(perms, orders, message): + """Every precondition is reported as ``ValueError``, never as a panic.""" with pytest.raises(ValueError, match=message): TranslationGroup.from_generators(4, perms, orders) +@pytest.mark.parametrize( + "ctor, args, message", + [ + (TranslationGroup.chain_1d, (0,), "n must be positive"), + (TranslationGroup.torus_2d, (0, 2), "lx must be positive"), + (TranslationGroup.torus_2d, (2, 0), "ly must be positive"), + (TranslationGroup.torus_3d, (2, 0, 2), "ly must be positive"), + (TranslationGroup.torus_3d, (2, 2, 0), "lz must be positive"), + (TranslationGroup.ladder, (0, 2), "l must be positive"), + (TranslationGroup.ladder, (2, 0), "n_legs must be positive"), + ], +) +def test_lattice_constructors_reject_empty_extents(ctor, args, message): + """Degenerate lattice extents raise ``ValueError`` rather than tripping + the core's ``assert!`` (which would surface as a ``PanicException``).""" + with pytest.raises(ValueError, match=message): + ctor(*args) + + def test_canonicalize_is_orbit_invariant(): g = TranslationGroup.chain_1d(4) shifts = ["IIXY", "IXYI", "XYII", "YIIX"] @@ -139,3 +166,16 @@ def test_symmetry_merge_rejects_qubit_count_mismatch(): p = psum(4, [("ZIII", 1.0)]) with pytest.raises(ValueError, match="4 qubits but the TranslationGroup acts on 3"): p.symmetry_merge(TranslationGroup.chain_1d(3)) + + +def test_lossy_pauli_sum_rejects_symmetry_merging(): + """`LossyPauliSum` inherits the merge wrappers but the Rust core expands + them only for non-loss variants, so both must fail with a clear + `NotImplementedError` rather than an `AttributeError` from inside the + wrapper.""" + lossy = LossyPauliSum(["ZIZI"], 4, [1.0]) + group = TranslationGroup.chain_1d(4) + with pytest.raises(NotImplementedError, match="not implemented for LossyPauliSum"): + lossy.symmetry_merge(group) + with pytest.raises(NotImplementedError, match="not implemented for LossyPauliSum"): + lossy.momentum_merge(LossyPauliSum(["ZIZI"], 4, [1.0]), group, [0])