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/basis.rs b/crates/ppvm-lindblad/src/basis.rs index 9f352d972..153780213 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,98 @@ impl LindbladSpec { } Ok(merged.into_iter().filter(|(_, c)| c.norm() > 0.0).collect()) } + + /// Phase-aware leakage: out-of-basis component of `L*(O_k)` where + /// `O_k` is the operator represented by `basis` (orbit reps) and + /// `coeffs` (complex coefficients in momentum `sector`). + /// + /// For each input rep `r` with coefficient `c_r`, and each output `q` + /// of `L*(r) = Σ_q v_q · q`: + /// 1. Canonicalize `q` → `(r_q, χ_k, |orbit_q|)` via + /// [`Sector::canonicalize_phase`]. + /// 2. If `r_q` NOT in `basis` and NOT in `protected`: + /// `merged[r_q] += χ_k · v_q · c_r · |orbit_r| / |orbit_q|`. + /// + /// The `|orbit_r| / |orbit_q|` factor is the convention conversion + /// documented on [`crate::mf_expm`]'s `build_orbit_rep_cols`, so the + /// admitted rates are comparable to the averaged-convention + /// coefficients the caller holds. + /// + /// Returns `(r_q, sum)` pairs for all candidates with nonzero sum. + /// + /// This is the orbit-rep counterpart of [`Self::leakage_with_prune`], + /// and caps the live candidate map the same way: to the *available + /// room* `room = max_basis − basis.len()` (the reps we could actually + /// add), applied during accumulation. A large `max_basis` + /// (room ≥ all candidates) disables the cap — the near-exact case. + pub fn leakage_orbit_rep( + &self, + basis: &[Word], + coeffs: &[Complex], + 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]; + // A rep that cannot carry the sector contributes + // nothing (its coefficient is identically zero). + let Some(orbit_in) = sector.orbit_size(r) else { + return Vec::new(); + }; + let terms = self.compute_action_terms(r, s1, s2, lm); + let mut out = Vec::with_capacity(terms.len()); + for (q, v) in terms.iter() { + let Some((r_q, phase, orbit_out)) = sector.canonicalize_phase(q) else { + continue; + }; + if !in_basis.contains(&r_q) && !protected_set.contains(&r_q) { + let rate = phase * *v * c_r * (orbit_in as f64 / orbit_out as f64); + out.push((r_q, rate)); + } + } + out + }, + ) + .collect(); + for v in local { + for (k, val) in v { + *merged.entry(k).or_insert(Complex::new(0.0, 0.0)) += val; + } + } + cap_map_to_room(&mut merged, room); + } + Ok(merged.into_iter().filter(|(_, c)| c.norm() > 0.0).collect()) + } } diff --git a/crates/ppvm-lindblad/src/lib.rs b/crates/ppvm-lindblad/src/lib.rs index da5414e66..57dff8aa2 100644 --- a/crates/ppvm-lindblad/src/lib.rs +++ b/crates/ppvm-lindblad/src/lib.rs @@ -38,8 +38,11 @@ mod basis; pub mod config; pub mod error; pub(crate) mod expm; +mod scalar; +pub mod sector; mod spec; mod step; +mod truncate; mod word; /// Matrix-free / quspin-expm-backed `exp(dt·L*)·b` engine. See module docs. @@ -48,6 +51,7 @@ pub(crate) mod mf_expm; pub use basis::build_basis_index; pub use config::PcStepConfig; pub use error::Error; +pub use sector::{Sector, canonicalize_basis_to_rep}; pub use spec::{JumpInput, LindbladSpec}; pub use step::PcStepTimings; pub use word::{MAX_QUBITS, Word, codes_from_word, parse_pauli_string, word_from_codes}; diff --git a/crates/ppvm-lindblad/src/mf_expm.rs b/crates/ppvm-lindblad/src/mf_expm.rs index f2e5b73f3..9437396f8 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,80 @@ 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 · |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 +/// 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)| { + // 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 Some((r_q, phase, orbit_out)) = sector.canonicalize_phase(q) else { + continue; + }; + if let Some(&row) = index.get(&r_q) { + let val = phase * *v * (orbit_in as f64 / orbit_out as f64); + 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 +296,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 +355,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 +368,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/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..b76ad2687 --- /dev/null +++ b/crates/ppvm-lindblad/src/sector.rs @@ -0,0 +1,106 @@ +// 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 · |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. +//! +//! ## 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`, 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) -> 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); + 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) + } +} + +/// 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..7a168a77e 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,92 @@ 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. + /// + /// Honours `cfg.num_threads` the same way [`Self::pc_step`] does. + pub fn pc_step_orbit_rep( + &self, + basis: &mut Vec, + coeffs: &mut Vec>, + dt: f64, + 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, + 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 4217b3e13..85526d37f 100644 --- a/crates/ppvm-lindblad/src/tests.rs +++ b/crates/ppvm-lindblad/src/tests.rs @@ -94,6 +94,147 @@ fn word_codec_roundtrip() { assert_eq!(out.as_slice(), &codes); } +/// 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'] { + let mut s = vec!['I'; n]; + s[j] = op; + s[nxt] = op; + h_terms.push((s.into_iter().collect(), 1.0)); + } + } + h_terms +} + +/// 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, 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 { + pc_step_complex_full(&spec, &mut bf, &mut cf, dt); + } + canonicalize_pauli_sum_complex(&mut bf, &mut cf, &group, k); + + // ----- 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); + let sector = Sector::new(&group, k); + for _ in 0..n_steps { + spec.pc_step_orbit_rep( + &mut br, + &mut cr, + dt, + &protected, + sector, + &PcStepConfig { + max_basis: 10_000_000, + ..Default::default() + }, + ) + .unwrap(); + } + + 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 {w} in orbit-rep but not in full-basis")); + 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}" + ); +} + +/// 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-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-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 30d6d4c01..b8ee315e1 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 //! @@ -66,9 +72,12 @@ 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}; +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..501cbceaa 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; @@ -50,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 @@ -99,6 +148,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(); @@ -122,24 +208,97 @@ 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); + projected.insert(rep, (sum, orbit_size)); } - basis.clear(); - coeffs.clear(); - basis.reserve(projected.len()); - coeffs.reserve(projected.len()); - for (w, c) in projected { - basis.push(w); - coeffs.push(c); + projected +} + +/// 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`]. +/// +/// 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 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 projected = project_onto_reps(&combined, group, k_modes); + re.data_mut().clear(); + im.data_mut().clear(); + for (word, (sum, _orbit_size)) in projected { + if sum.re != 0.0 { + *re += (word, sum.re); + } + 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 64c1deec0..66c261d3f 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; @@ -135,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); @@ -452,6 +493,159 @@ 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` 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; + use crate::prelude::*; + + type Cfg = ByteFxHashF64<1>; + let n = 4usize; + let group = TranslationGroup::chain_1d(n); + + 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'; + 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); + } + } + // 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); + + 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() { @@ -472,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/interface.rs b/crates/ppvm-python-native/src/interface.rs index 44a47b83d..6657b5451 100644 --- a/crates/ppvm-python-native/src/interface.rs +++ b/crates/ppvm-python-native/src/interface.rs @@ -48,6 +48,104 @@ 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.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, 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 + /// 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, + 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()), + ("other", other.inner.n_qubits()), + ] { + if n != n_q { + return Err(pyo3::exceptions::PyValueError::new_err(format!( + "{label} PauliSum has {n} qubits but the \ + TranslationGroup acts on {n_q}", + ))); + } + } + 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(), + ))); + } + ppvm_pauli_sum::symmetry::momentum_merge_pauli_sum_pair( + &mut self.inner, + &mut other.inner, + group.core(), + &momentum, + ); + Ok(()) + } + } + }; +} + macro_rules! create_strategy { (false, $min_abs_coeff:ident, $max_pauli_weight:ident, $_max_loss_weight:ident) => { CombinedStrategy( @@ -403,6 +501,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..e4f9228c8 100644 --- a/crates/ppvm-python-native/src/lib.rs +++ b/crates/ppvm-python-native/src/lib.rs @@ -15,7 +15,9 @@ 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; pub(crate) fn flat_pairs(targets: &[usize]) -> PyResult> { if !targets.len().is_multiple_of(2) { @@ -310,4 +312,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..477e1e834 100644 --- a/crates/ppvm-python-native/src/lindblad.rs +++ b/crates/ppvm-python-native/src/lindblad.rs @@ -13,18 +13,19 @@ use std::collections::HashMap; use num::Complex; -use numpy::{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>); +type PyPauliMapComplex<'py> = (Bound<'py, PyArray2>, Bound<'py, PyArray1>); type PyCoo<'py> = ( Bound<'py, PyArray1>, Bound<'py, PyArray1>, 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()) } @@ -43,29 +44,9 @@ 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_group_qubits, check_momentum_len, decode_basis, encode_basis, +}; /// Pack `Vec<(Word, f64)>` into the standard PyO3 return shape. fn pack_pauli_map<'py>( @@ -73,17 +54,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))) } @@ -175,13 +147,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)? @@ -238,13 +204,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 { @@ -300,13 +260,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 { @@ -341,6 +295,106 @@ 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, + num_threads = 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, + num_threads: Option, + ) -> PyResult> { + use num::Complex; + use ppvm_lindblad::{Sector, canonicalize_basis_to_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()?; + 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)) + .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()?; + 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, + &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, + }, + ) + .map_err(map_err)?; + + let out_coeffs: Vec = coeffs_vec + .iter() + .map(|c| Complex64::new(c.re, c.im)) + .collect(); + let basis_arr = encode_basis(py, &basis_words, n_q)?; + 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/pauli_arr.rs b/crates/ppvm-python-native/src/pauli_arr.rs new file mode 100644 index 000000000..45580d638 --- /dev/null +++ b/crates/ppvm-python-native/src/pauli_arr.rs @@ -0,0 +1,107 @@ +// 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 [`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 { + 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 new file mode 100644 index 000000000..cd3dfa32b --- /dev/null +++ b/crates/ppvm-python-native/src/symmetry.rs @@ -0,0 +1,282 @@ +// 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, 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>); + +/// 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 + } +} + +/// 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) -> PyResult { + check_lattice(&[("n", n)])?; + Ok(Self { + inner: core_sym::TranslationGroup::chain_1d(n), + }) + } + + #[staticmethod] + 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) -> 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) -> 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 { + 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. + #[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/|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 +/// [`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(); + 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()?; + 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)) + .collect(); + + core_sym::canonicalize_pauli_sum_complex( + &mut basis_words, + &mut coeffs_vec, + &group.inner, + k_slice, + ); + + let out_coeffs: Vec = coeffs_vec + .iter() + .map(|c| Complex64::new(c.re, c.im)) + .collect(); + let basis_arr = encode_basis(py, &basis_words, n_q)?; + 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(); + 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()?; + 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)) + .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(); + check_group_width(&basis_view, n_q)?; + let coeffs_slice = coeffs.as_slice()?; + check_coeffs_len(coeffs_slice.len(), basis_view.shape()[0])?; + + 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); + + 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/_core.pyi b/ppvm-python/src/ppvm/_core.pyi index 4b2cfac7f..e9877b7f8 100644 --- a/ppvm-python/src/ppvm/_core.pyi +++ b/ppvm-python/src/ppvm/_core.pyi @@ -63,6 +63,17 @@ class _PauliSumBase: def weights(self) -> list[tuple[str, int]]: ... def current_max_weight(self) -> int: ... +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: _PauliSumNoLossBase, + group: TranslationGroup, + momentum: list[int], + ) -> None: ... + class _PauliSumLossBase(_PauliSumBase): def loss_channel(self, addr0: int, p: float, truncate: bool = True) -> None: ... def correlated_loss_channel( @@ -70,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): ... @@ -392,4 +403,57 @@ 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, + num_threads: int | 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: ... diff --git a/ppvm-python/src/ppvm/lindblad.py b/ppvm-python/src/ppvm/lindblad.py index a31ed4e5f..a71dce745 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 @@ -41,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} @@ -290,6 +293,76 @@ 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: _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. + + 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. + + 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 + 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). + + ``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: + 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), + None if num_threads is None else int(num_threads), + ) + 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..6d1922915 100644 --- a/ppvm-python/src/ppvm/paulisum.py +++ b/ppvm-python/src/ppvm/paulisum.py @@ -386,6 +386,60 @@ def trace(self, pattern: str) -> float: """ return self._interface.trace(pattern) + 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) + 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.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: _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 + 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. 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 + error as the ``k=0`` merge. + + Args: + other: the PauliSum holding the imaginary component (modified in place). + 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). + """ + 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. @@ -497,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/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 new file mode 100644 index 000000000..6bb8a81d2 --- /dev/null +++ b/ppvm-python/test/lindblad/test_pc_step_orbit_rep.py @@ -0,0 +1,320 @@ +# 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, 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): + """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): + """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_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.""" + 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} + + +@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) + 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)) + 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(): + 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 new file mode 100644 index 000000000..5f8c02a70 --- /dev/null +++ b/ppvm-python/test/test_momentum_merge.py @@ -0,0 +1,245 @@ +# 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, 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 + + +@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 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) + 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) + # 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) + 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 + 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 diff --git a/ppvm-python/test/test_symmetry_arrays.py b/ppvm-python/test/test_symmetry_arrays.py new file mode 100644 index 000000000..7630e6cc0 --- /dev/null +++ b/ppvm-python/test/test_symmetry_arrays.py @@ -0,0 +1,242 @@ +# 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``, 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 + 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 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 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 + 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..f91536534 --- /dev/null +++ b/ppvm-python/test/test_symmetry_merge.py @@ -0,0 +1,181 @@ +# 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 LossyPauliSum, PauliSum, 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"), + # 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"] + 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)) + + +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])