Adaptive Pauli-Lindbladian shim + adaptive-evolution demo

Cargo.toml

@@ -14,6 +14,7 @@ ppvm-sym = { version = "0.1.0", path = "crates/ppvm-sym" }
 members = [
     "crates/ppvm-runtime",
     "crates/ppvm-sym",
+    "crates/ppvm-lindblad",
     "crates/ppvm-python-native",
     "crates/ppvm-tableau",
     "crates/ppvm-stim", "crates/stim-parser",

crates/ppvm-lindblad/Cargo.toml

@@ -0,0 +1,24 @@
+[package]
+name = "ppvm-lindblad"
+version = "0.1.0"
+edition = "2024"
+description = "Direct Heisenberg-picture Lindbladian evolution on an adaptive Pauli-string basis."
+
+[dependencies]
+fxhash = "0.2.1"
+ndarray = "0.17"
+num = "0.4.3"
+ppvm-runtime = { version = "0.1.0", path = "../ppvm-runtime" }
+rayon = "1.11"
+sprs = { version = "0.11", default-features = false, features = ["multi_thread"] }
+# Matrix-exponential action (Al-Mohy & Higham). NOTE: QuSpin-rust currently has
+# no LICENSE file; this pin is a placeholder pending a permissive license being
+# added upstream (the authors are co-authors of the ppvm paper).
+quspin-expm = { git = "https://github.com/QuSpin/QuSpin-rust", rev = "6f4a1581e14d5258c9aeb684b4aec964bdff4e34" }
+# `QuSpinError` (the error type returned by the `LinearOperator` trait methods
+# 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 = "6f4a1581e14d5258c9aeb684b4aec964bdff4e34" }
+
+[dev-dependencies]
+approx = "0.5.1"

crates/ppvm-lindblad/src/config.rs

@@ -0,0 +1,27 @@
+// SPDX-FileCopyrightText: 2026 The PPVM Authors
+// SPDX-License-Identifier: Apache-2.0
+
+//! Configuration objects for the predictor-corrector stepper.
+
+/// Policy options for a single predictor-corrector step
+/// ([`crate::LindbladSpec::pc_step`] / [`crate::LindbladSpec::pc_step_timed`]).
+///
+/// These are the per-run *tuning knobs*, kept separate from the per-call data
+/// (`basis`, `coeffs`, `dt`, `protected`) so the step functions don't grow an
+/// unwieldy positional-argument list.
+#[derive(Debug, Clone, Copy)]
+pub struct PcStepConfig {
+    /// Leakage-rate threshold for basis growth: a leakage Pauli is appended to
+    /// the basis only if its `|coeff|` exceeds `tau_add`. Larger values keep the
+    /// basis smaller (faster, less accurate); it also drives the streaming
+    /// Cauchy-Schwarz prune in [`crate::LindbladSpec::leakage_with_prune`].
+    pub tau_add: f64,
+    /// Magnitude prune applied after the corrector: basis entries whose
+    /// `|coeff|` is below `drop_tol` are discarded (protected words are always
+    /// kept). `drop_tol <= 0.0` disables pruning.
+    pub drop_tol: f64,
+    /// When `Some(n)`, run the entire step inside a freshly built rayon thread
+    /// pool of `n` threads (useful for benchmarking parallel scaling). When
+    /// `None`, the global rayon pool is used.
+    pub num_threads: Option<usize>,
+}

crates/ppvm-lindblad/src/error.rs

@@ -0,0 +1,76 @@
+// SPDX-FileCopyrightText: 2026 The PPVM Authors
+// SPDX-License-Identifier: Apache-2.0
+
+//! Error type for [`crate::LindbladSpec`] construction and stepping.
+
+use crate::MAX_QUBITS;
+use std::fmt;
+
+/// Errors raised when constructing a [`crate::LindbladSpec`].
+#[derive(Debug, Clone)]
+pub enum Error {
+    TooManyQubits {
+        got: usize,
+    },
+    LengthMismatch {
+        what: &'static str,
+        a: usize,
+        b: usize,
+    },
+    InvalidPauliCode {
+        code: u8,
+    },
+    InvalidPauliChar {
+        c: char,
+    },
+    WrongLength {
+        expected: usize,
+        got: usize,
+    },
+    NegativeRate {
+        index: usize,
+        rate: f64,
+    },
+    EmptyLincomb {
+        index: usize,
+    },
+    Internal(String),
+}
+
+impl fmt::Display for Error {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            Error::TooManyQubits { got } => {
+                write!(
+                    f,
+                    "LindbladSpec supports n_qubits ≤ {MAX_QUBITS}; got {got}"
+                )
+            }
+            Error::LengthMismatch { what, a, b } => {
+                write!(f, "{what}: expected matching lengths, got {a} and {b}")
+            }
+            Error::InvalidPauliCode { code } => write!(
+                f,
+                "Pauli code must be 0 (I), 1 (X), 2 (Z), or 3 (Y); got {code}"
+            ),
+            Error::InvalidPauliChar { c } => {
+                write!(f, "invalid Pauli character '{c}'; expected I, X, Y, or Z")
+            }
+            Error::WrongLength { expected, got } => {
+                write!(f, "Pauli string has length {got} but n_qubits = {expected}")
+            }
+            Error::NegativeRate { index, rate } => {
+                write!(f, "jump rate must be non-negative; got γ_{index} = {rate}")
+            }
+            Error::EmptyLincomb { index } => {
+                write!(
+                    f,
+                    "jump {index}: lincomb must contain at least one Pauli term"
+                )
+            }
+            Error::Internal(msg) => write!(f, "internal error: {msg}"),
+        }
+    }
+}
+
+impl std::error::Error for Error {}

crates/ppvm-lindblad/src/expm.rs

@@ -0,0 +1,99 @@
+// SPDX-FileCopyrightText: 2026 The PPVM Authors
+// SPDX-License-Identifier: Apache-2.0
+
+//! Support pieces for the `quspin-expm`-backed `exp(t·A)·b` engine.
+//!
+//! The hand-rolled scaling-and-squaring engine that once lived here has been
+//! retired in favour of the external `quspin-expm` crate (driven from
+//! [`crate::mf_expm`]). What remains is the shared scaffolding the new path
+//! still relies on: the `(m, s)` selection table [`THETA`] / [`select_ms`]
+//! from Al-Mohy & Higham (2011), used to pick the Taylor partition handed to
+//! `quspin-expm`'s `from_parts`.
+
+/// `θ_m` table from Al-Mohy & Higham (2011), Table A.3, for double
+/// precision (unit roundoff `u = 2^{-53}`).
+///
+/// `θ_m` bounds `‖A‖₁` such that the degree-`m` Taylor polynomial
+/// approximates `exp(A)` to within `u`. We pick `(m, s)` with
+/// `s ≥ ⌈‖tA‖₁ / θ_m⌉` and minimise `m·s` (total SpMV count).
+pub(crate) const THETA: &[(u32, f64)] = &[
+    (1, 2.29e-16),
+    (2, 2.58e-8),
+    (3, 1.39e-5),
+    (4, 3.40e-4),
+    (5, 2.40e-3),
+    (6, 9.07e-3),
+    (7, 2.38e-2),
+    (8, 5.00e-2),
+    (9, 8.96e-2),
+    (10, 1.44e-1),
+    (11, 2.14e-1),
+    (12, 3.00e-1),
+    (13, 4.00e-1),
+    (14, 5.14e-1),
+    (15, 6.41e-1),
+    (16, 7.81e-1),
+    (17, 9.31e-1),
+    (18, 1.09),
+    (19, 1.26),
+    (20, 1.44),
+    (21, 1.62),
+    (22, 1.82),
+    (23, 2.01),
+    (24, 2.22),
+    (25, 2.43),
+    (26, 2.64),
+    (27, 2.86),
+    (28, 3.08),
+    (29, 3.31),
+    (30, 3.54),
+];
+
+/// Pick `(m, s)` minimising `s·m` subject to `s ≥ ⌈t_norm / θ_m⌉, s ≥ 1`.
+/// Restricted to `m ≤ 30`; for larger norms `s` simply grows linearly.
+/// When `max_m` is set, only entries with `m ≤ max_m` are considered.
+pub(crate) fn select_ms(t_norm: f64, max_m: Option<u32>) -> (u32, u32) {
+    if t_norm <= 0.0 {
+        return (1, 1);
+    }
+    let mut best_m = 1u32;
+    let mut best_s = 1u32;
+    let mut best_cost = u64::MAX;
+    for &(m, theta) in THETA {
+        if let Some(cap) = max_m {
+            if m > cap {
+                continue;
+            }
+        }
+        let s_f = (t_norm / theta).ceil();
+        let s = if s_f >= 1.0 { s_f as u32 } else { 1 };
+        let cost = (m as u64) * (s as u64);
+        if cost < best_cost {
+            best_cost = cost;
+            best_m = m;
+            best_s = s;
+        }
+    }
+    (best_m, best_s)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn ms_selection_sane() {
+        // tiny norm → small m, s = 1
+        let (m, s) = select_ms(1e-9, None);
+        assert!(m <= 5, "expected small m for tiny norm, got m={m}");
+        assert_eq!(s, 1);
+
+        // moderate norm → m·s should be ~10-50
+        let (m, s) = select_ms(1.0, None);
+        assert!((m * s) <= 50, "moderate norm cost too high: m={m} s={s}");
+
+        // large norm → s grows
+        let (_m, s) = select_ms(100.0, None);
+        assert!(s >= 20, "large norm should require many steps, got s={s}");
+    }
+}