wai_quantum/

quantum_cal.rs

//! Deterministic autonomous qubit CALIBRATION — `wai.quantum.cal.rabi`
//! (extensions/quantum-ops § Calibration engines).
//!
//! **Autonomous gate tune-up** — sweep a control parameter, measure the qubit, fit
//! the response, set the optimum, no human in the loop — is a set of *algorithms
//! that run on measured data*, commonly developed against **device simulators**.
//! This module is that algorithm — the canonical **Rabi π-pulse tune-up** — done
//! the WAI way: a deterministic fixed-point device simulator + an autonomous
//! calibrator, **byte-identical on every machine** (a portable tune-up hash), so
//! a calibration run reproduces exactly and drops straight into a signed,
//! joule-metered [`crate::quantum_ops::CalibrationReceipt`].
//!
//! Two honest boundaries: (1) it tunes a **simulated** device — we have no
//! cryostat, exactly the surface a digital-twin demo runs on; (2) the physics
//! model is a clean Rabi model, not a full transmon Hamiltonian. What is *not*
//! approximate is the reproducibility: no float anywhere in the run, so the
//! calibrated amplitude and the outcome hash are the same bytes everywhere — the
//! property a float run does not have, on top of which the receipt
//! makes every run auditable.

/// Fixed-point fractional bits. All amplitudes, probabilities and angles are i64
/// at `2^CAL_FRAC`. Probabilities live in `[0, CAL_ONE]`.
pub const CAL_FRAC: u32 = 20;
/// `1.0` in fixed-point.
pub const CAL_ONE: i64 = 1 << CAL_FRAC;
/// `π · 2^CAL_FRAC` (integer constant — no float at runtime).
const PI_FX: i64 = 3_294_199; // round(π · 2^20)

#[inline]
fn fmul(a: i64, b: i64) -> i64 {
    ((a as i128 * b as i128) >> CAL_FRAC) as i64
}

/// Fixed-point sine via the Bhaskara-I rational approximation (max error ≈ 0.16%,
/// pure integer). Reduced modulo `π`: callers square the result, and `sin²` has
/// period `π`, so the sign dropped by the reduction is irrelevant. Deterministic
/// by construction — the same bytes on every machine.
pub fn sin2_fx(theta: i64) -> i64 {
    let pi = PI_FX as i128;
    let t = (theta as i128).rem_euclid(pi); // [0, π)
    let p = t * (pi - t); // t(π−t) at 2^(2·CAL_FRAC)
    let num = 16 * p;
    let den = 5 * pi * pi - 4 * p;
    let s = ((num << CAL_FRAC) / den) as i64; // sin(t) at 2^CAL_FRAC, in [0, ~CAL_ONE]
    fmul(s, s) // sin²(t) at 2^CAL_FRAC
}

// ---------------------------------------------------------------------------
// Deterministic device simulator
// ---------------------------------------------------------------------------

fn splitmix64(state: &mut u64) -> u64 {
    *state = state.wrapping_add(0x9E37_79B9_7F4A_7C15);
    let mut z = *state;
    z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
    z ^ (z >> 31)
}

/// A single driven qubit under a Rabi model. The calibrator does not see
/// `a_pi_fx` — it is the hidden truth the tune-up must discover.
#[derive(Clone, Copy, Debug)]
pub struct Device {
    /// The true π-pulse amplitude (fixed-point) — hidden from the calibrator.
    pub a_pi_fx: i64,
    /// A per-device stream seed, so two devices with the same `a_pi` still shot-
    /// noise independently but each reproducibly.
    pub seed: u64,
}

impl Device {
    /// Exact excited-state probability at drive amplitude `a`:
    /// `P(a) = sin²(π·a / (2·a_π))`. Pure integer.
    pub fn excited_prob(&self, a_fx: i64) -> i64 {
        // θ = π·a / (2·a_π), fixed-point
        let theta = ((PI_FX as i128 * a_fx as i128) / (2 * self.a_pi_fx as i128)) as i64;
        sin2_fx(theta)
    }

    /// One "experiment": `shots` projective measurements at amplitude `a`, returning
    /// the measured excited fraction (fixed-point). Deterministic Bernoulli sampling
    /// via `splitmix64` keyed by `(seed, a)`, so the shot noise is real but portable.
    pub fn measure(&self, a_fx: i64, shots: u32) -> i64 {
        let p = self.excited_prob(a_fx);
        let mut st = self
            .seed
            .wrapping_mul(0x1000_0001)
            .wrapping_add(a_fx as u64);
        let mut hits: u64 = 0;
        for _ in 0..shots {
            let r = (splitmix64(&mut st) >> (64 - CAL_FRAC)) as i64; // [0, CAL_ONE)
            if r < p {
                hits += 1;
            }
        }
        ((hits as i128 * CAL_ONE as i128) / shots.max(1) as i128) as i64
    }
}

// ---------------------------------------------------------------------------
// Autonomous Rabi tune-up
// ---------------------------------------------------------------------------

/// Sweep + refinement schedule for [`calibrate_rabi`].
#[derive(Clone, Copy, Debug)]
pub struct TuneConfig {
    /// Upper bound of the amplitude sweep (fixed-point). Choose ≳ 1.5× the
    /// expected π-amplitude so the first Rabi lobe is contained.
    pub a_max_fx: i64,
    /// Coarse-sweep points across `[0, a_max]`.
    pub n_coarse: u32,
    /// Fine-sweep points around the coarse peak (± one coarse step).
    pub n_fine: u32,
    /// Shots per measurement.
    pub shots: u32,
}

impl Default for TuneConfig {
    fn default() -> Self {
        TuneConfig { a_max_fx: 2 * CAL_ONE, n_coarse: 41, n_fine: 41, shots: 2000 }
    }
}

/// The result of an autonomous tune-up. Deterministic given `(device, cfg)`.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct TuneResult {
    /// The calibrated π-pulse amplitude (fixed-point).
    pub a_pi_est_fx: i64,
    /// Measured excited fraction at the calibrated amplitude — the tune-up's own
    /// evidence of a good π rotation (`≈ CAL_ONE` is a clean π-pulse).
    pub peak_prob_fx: i64,
    /// Total measurements performed (the deterministic "work" of the run).
    pub n_measurements: u32,
    /// The full sweep: `(amplitude, measured_prob)` at `2^CAL_FRAC`, coarse then fine.
    pub sweep: Vec<(i64, i64)>,
}

impl TuneResult {
    /// BLAKE3 over the canonical outcome — a portable tune-up hash. Two machines
    /// running the same `(device, cfg)` get this same hash.
    pub fn hash(&self) -> [u8; 32] {
        let mut h = blake3::Hasher::new();
        h.update(b"wai:quantum-cal-rabi\x01");
        h.update(&self.a_pi_est_fx.to_le_bytes());
        h.update(&self.peak_prob_fx.to_le_bytes());
        h.update(&(self.n_measurements as u64).to_le_bytes());
        for (a, p) in &self.sweep {
            h.update(&a.to_le_bytes());
            h.update(&p.to_le_bytes());
        }
        *h.finalize().as_bytes()
    }

    /// Vendor-neutral artifacts for [`seal_artifacts`] — uniform with the other
    /// routines' `artifacts()`.
    pub fn artifacts(&self) -> CalArtifacts {
        let mut ev = Vec::with_capacity(self.sweep.len() * 16 + 16);
        ev.extend_from_slice(&self.a_pi_est_fx.to_le_bytes());
        ev.extend_from_slice(&self.peak_prob_fx.to_le_bytes());
        for (a, p) in &self.sweep {
            ev.extend_from_slice(&a.to_le_bytes());
            ev.extend_from_slice(&p.to_le_bytes());
        }
        CalArtifacts {
            target: "1q_gate_amplitude".into(),
            config_bytes: self.a_pi_est_fx.to_le_bytes().to_vec(),
            evidence_kind: "rabi_amplitude_sweep".into(),
            evidence_bytes: ev,
            summary: format!("pi_amp={} peak_prob={} (2^{CAL_FRAC})", self.a_pi_est_fx, self.peak_prob_fx),
        }
    }
}

/// The autonomous calibrator: sweep the drive amplitude, then **fit the whole
/// Rabi curve to the model** `P(a) = sin²(π·a/2a_π)` — searching the candidate
/// π-amplitude that minimizes total error, coarse then fine. Whole-curve fitting
/// (what real tune-ups do) is robust where an argmax fails: near a *flat* maximum
/// shot noise moves the peak sample, but the model fit uses every point. Pure
/// integer — byte-identical on every machine.
pub fn calibrate_rabi(device: &Device, cfg: &TuneConfig) -> TuneResult {
    let n = cfg.n_coarse.max(4);
    let mut sweep = Vec::with_capacity(n as usize);
    for i in 0..n {
        let a = ((cfg.a_max_fx as i128 * i as i128) / (n - 1) as i128) as i64;
        sweep.push((a, device.measure(a, cfg.shots)));
    }

    // Model-fit a_π over a candidate grid: minimize Σ|measured − sin²(π·a/2a_π)|.
    let fit = |lo: i64, hi: i64, steps: u32| -> i64 {
        let steps = steps.max(2);
        let mut best = lo.max(1);
        let mut best_err = i128::MAX;
        for g in 1..=steps {
            let a_pi = lo + ((hi - lo) as i128 * g as i128 / steps as i128) as i64;
            if a_pi <= 0 {
                continue;
            }
            let mut err = 0i128;
            for &(a, p) in &sweep {
                let theta = ((PI_FX as i128 * a as i128) / (2 * a_pi as i128)) as i64;
                err += (p - sin2_fx(theta)).unsigned_abs() as i128;
            }
            if err < best_err {
                best_err = err;
                best = a_pi;
            }
        }
        best
    };
    let coarse_step = cfg.a_max_fx / (n as i64 - 1);
    let a_coarse = fit(coarse_step, cfg.a_max_fx, cfg.n_coarse.max(24));
    let a_pi = fit((a_coarse - coarse_step).max(1), a_coarse + coarse_step, cfg.n_fine.max(24));

    // A verification pulse at the calibrated amplitude — the π-rotation quality.
    let peak = device.measure(a_pi, cfg.shots);
    TuneResult { a_pi_est_fx: a_pi, peak_prob_fx: peak, n_measurements: n + 1, sweep }
}

// ---------------------------------------------------------------------------
// Calibration → signed receipt
// ---------------------------------------------------------------------------

#[cfg(feature = "quantum_ops")]
mod sealed {
    use super::*;
    use crate::quantum_ops::{content_hash, CalibrationReceipt, Evidence, GrantRef};
    use ed25519_dalek::SigningKey;

    /// Serialize the sweep to canonical bytes (the "raw dataset" the receipt
    /// content-addresses — the analog of a vendor HDF5 tuning trace).
    pub fn sweep_bytes(r: &TuneResult) -> Vec<u8> {
        let mut o = Vec::with_capacity(r.sweep.len() * 16 + 16);
        o.extend_from_slice(&r.a_pi_est_fx.to_le_bytes());
        o.extend_from_slice(&r.peak_prob_fx.to_le_bytes());
        for (a, p) in &r.sweep {
            o.extend_from_slice(&a.to_le_bytes());
            o.extend_from_slice(&p.to_le_bytes());
        }
        o
    }

    /// Run an autonomous Rabi tune-up and seal a `wai.quantum.calibration`
    /// receipt over it: the config hash pins the calibrated amplitude, the
    /// evidence content-addresses the full sweep, the joules are the run's
    /// measured energy, and the grant authorizes it. The tune-up is byte-exact,
    /// so a verifier re-running `(device, cfg)` reproduces the sealed evidence.
    #[allow(clippy::too_many_arguments)]
    pub fn calibrate_and_seal(
        signer: &SigningKey,
        signer_id: impl Into<String>,
        device_id: impl Into<String>,
        target: impl Into<String>,
        device: &Device,
        cfg: &TuneConfig,
        joules_micro: u64,
        grant: GrantRef,
        parent: Option<[u8; 32]>,
    ) -> (TuneResult, CalibrationReceipt) {
        let r = calibrate_rabi(device, cfg);
        // config = the tuned control parameter (the calibrated π-amplitude).
        let config_hash = content_hash(&r.a_pi_est_fx.to_le_bytes());
        // evidence = the sweep dataset + the achieved π-fidelity proxy.
        let ev = vec![Evidence {
            kind: "rabi_amplitude_sweep".into(),
            blob_hash: content_hash(&sweep_bytes(&r)),
            summary: format!(
                "pi_amp={} peak_prob={} (fixed-point 2^{})",
                r.a_pi_est_fx, r.peak_prob_fx, CAL_FRAC
            ),
        }];
        let receipt = CalibrationReceipt::seal(
            signer, signer_id, device_id, target, config_hash, ev, joules_micro, grant,
            parent,
        );
        (r, receipt)
    }
}

#[cfg(feature = "quantum_ops")]
pub use sealed::{calibrate_and_seal, sweep_bytes};

// ===========================================================================
// Shared fixed-point machinery for the other routines
// ===========================================================================

const TWO_PI_FX: i64 = 6_588_397; // round(2π · 2^20)
const HALF_PI_FX: i64 = 1_647_099; // round(π/2 · 2^20)
const LN2_FX: i64 = 726_817; // round(ln2 · 2^20)

/// Signed sine over a full period, `[-CAL_ONE, CAL_ONE]`, via Bhaskara on each half.
pub fn sin_fx(theta: i64) -> i64 {
    let two_pi = TWO_PI_FX as i128;
    let mut t = (theta as i128).rem_euclid(two_pi); // [0, 2π)
    let pi = PI_FX as i128;
    let sign = if t >= pi {
        t -= pi;
        -1i64
    } else {
        1
    };
    let p = t * (pi - t);
    let num = 16 * p;
    let den = 5 * pi * pi - 4 * p;
    let s = ((num << CAL_FRAC) / den) as i64;
    sign * s
}

/// Signed cosine, `cos(x) = sin(x + π/2)`.
pub fn cos_fx(theta: i64) -> i64 {
    sin_fx(theta + HALF_PI_FX)
}

/// `exp(-x)` for `x ≥ 0`, fixed-point. Splits `x = k·ln2 + r`, a cubic on the
/// small remainder, and `2^-k` by shift. Pure integer, deterministic.
pub fn exp_neg_fx(x_fx: i64) -> i64 {
    if x_fx <= 0 {
        return CAL_ONE;
    }
    let k = x_fx / LN2_FX;
    let r = x_fx - k * LN2_FX; // [0, ln2)
    let r2 = fmul(r, r);
    let r3 = fmul(r2, r);
    let er = CAL_ONE - r + r2 / 2 - r3 / 6; // exp(-r)
    if k >= 62 {
        return 0;
    }
    (er >> k).max(0)
}

fn ratio_fx(a: i64, b: i64) -> i64 {
    if b == 0 {
        return 0;
    }
    (((a as i128) << CAL_FRAC) / b as i128) as i64
}

/// The vendor-neutral artifacts of any calibration run — what a receipt binds:
/// the target, the calibrated parameter (`config`), and the evidence dataset.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct CalArtifacts {
    pub target: String,
    pub config_bytes: Vec<u8>,
    pub evidence_kind: String,
    pub evidence_bytes: Vec<u8>,
    pub summary: String,
}

// ===========================================================================
// Ramsey / T2* — dephasing + detuning
// ===========================================================================

/// A driven qubit under a Ramsey model: free-induction with an angular detuning
/// `omega` and a dephasing time `t2`. The calibrator sees neither.
#[derive(Clone, Copy, Debug)]
pub struct RamseyDevice {
    /// Angular detuning ω (rad per unit time, fixed-point) — hidden.
    pub omega_fx: i64,
    /// Dephasing time T2 (fixed-point) — hidden.
    pub t2_fx: i64,
    pub seed: u64,
}

impl RamseyDevice {
    /// `P(τ) = ½(1 + cos(ω·τ)·exp(-τ/T2))`.
    pub fn prob(&self, tau_fx: i64) -> i64 {
        let c = cos_fx(fmul(self.omega_fx, tau_fx));
        let env = exp_neg_fx(ratio_fx(tau_fx, self.t2_fx));
        (CAL_ONE + fmul(c, env)) / 2
    }
    pub fn measure(&self, tau_fx: i64, shots: u32) -> i64 {
        let p = self.prob(tau_fx);
        let mut st = self.seed.wrapping_mul(0x2545_F491).wrapping_add(tau_fx as u64);
        let mut hits = 0u64;
        for _ in 0..shots {
            let r = (splitmix64(&mut st) >> (64 - CAL_FRAC)) as i64;
            if r < p {
                hits += 1;
            }
        }
        ((hits as i128 * CAL_ONE as i128) / shots.max(1) as i128) as i64
    }
}

/// Ramsey autonomous fit — coarse search for the detuning ω (correlation) then
/// the dephasing T2 (least model error), both over deterministic parameter grids.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct RamseyResult {
    pub omega_est_fx: i64,
    pub t2_est_fx: i64,
    pub sweep: Vec<(i64, i64)>,
    pub n_measurements: u32,
}

impl RamseyResult {
    pub fn hash(&self) -> [u8; 32] {
        let mut h = blake3::Hasher::new();
        h.update(b"wai:quantum-cal-ramsey\x01");
        h.update(&self.omega_est_fx.to_le_bytes());
        h.update(&self.t2_est_fx.to_le_bytes());
        for (a, b) in &self.sweep {
            h.update(&a.to_le_bytes());
            h.update(&b.to_le_bytes());
        }
        *h.finalize().as_bytes()
    }
    pub fn artifacts(&self) -> CalArtifacts {
        let mut cfg = Vec::new();
        cfg.extend_from_slice(&self.omega_est_fx.to_le_bytes());
        cfg.extend_from_slice(&self.t2_est_fx.to_le_bytes());
        let mut ev = Vec::new();
        for (a, b) in &self.sweep {
            ev.extend_from_slice(&a.to_le_bytes());
            ev.extend_from_slice(&b.to_le_bytes());
        }
        CalArtifacts {
            target: "dephasing".into(),
            config_bytes: cfg,
            evidence_kind: "ramsey_fringe".into(),
            evidence_bytes: ev,
            summary: format!("omega={} t2={} (2^{CAL_FRAC})", self.omega_est_fx, self.t2_est_fx),
        }
    }
}

/// `tau_max`, points, shots; `omega_max`/`t2_max` bound the fit search grids.
#[derive(Clone, Copy, Debug)]
pub struct RamseyConfig {
    pub tau_max_fx: i64,
    pub n_points: u32,
    pub shots: u32,
    pub omega_max_fx: i64,
    pub t2_max_fx: i64,
    pub n_grid: u32,
}

pub fn calibrate_ramsey(dev: &RamseyDevice, cfg: &RamseyConfig) -> RamseyResult {
    let n = cfg.n_points.max(4);
    let mut sweep = Vec::with_capacity(n as usize);
    for i in 0..n {
        let tau = ((cfg.tau_max_fx as i128 * i as i128) / (n - 1) as i128) as i64;
        sweep.push((tau, dev.measure(tau, cfg.shots)));
    }
    // detuning: maximize |Σ (P−½)·cos(ω·τ)| over candidate ω
    let ng = cfg.n_grid.max(4);
    let mut best_omega = 0i64;
    let mut best_score = -1i128;
    for g in 1..ng {
        let omega = ((cfg.omega_max_fx as i128 * g as i128) / ng as i128) as i64;
        let mut acc = 0i128;
        for &(tau, p) in &sweep {
            acc += ((p - CAL_ONE / 2) as i128 * cos_fx(fmul(omega, tau)) as i128) >> CAL_FRAC;
        }
        let s = acc.abs();
        if s > best_score {
            best_score = s;
            best_omega = omega;
        }
    }
    // T2: minimize Σ |P − model(ω_est, T2)| over candidate T2
    let mut best_t2 = cfg.t2_max_fx;
    let mut best_err = i128::MAX;
    for g in 1..=ng {
        let t2 = ((cfg.t2_max_fx as i128 * g as i128) / ng as i128) as i64;
        let probe = RamseyDevice { omega_fx: best_omega, t2_fx: t2, seed: 0 };
        let mut err = 0i128;
        for &(tau, p) in &sweep {
            err += (p - probe.prob(tau)).unsigned_abs() as i128;
        }
        if err < best_err {
            best_err = err;
            best_t2 = t2;
        }
    }
    RamseyResult { omega_est_fx: best_omega, t2_est_fx: best_t2, sweep, n_measurements: n }
}

// ===========================================================================
// Randomized benchmarking — average gate fidelity
// ===========================================================================

/// A gate set with a per-gate depolarizing error `epsilon` — hidden.
#[derive(Clone, Copy, Debug)]
pub struct RbDevice {
    pub epsilon_fx: i64,
    pub seed: u64,
}

impl RbDevice {
    /// Survival probability of a length-`m` Clifford sequence:
    /// `P(m) = ½ + ½·p^m`, `p = 1 − 2ε`.
    pub fn survival(&self, m: u32) -> i64 {
        let p = CAL_ONE - 2 * self.epsilon_fx;
        let mut pm = CAL_ONE;
        for _ in 0..m {
            pm = fmul(pm, p);
        }
        CAL_ONE / 2 + fmul(CAL_ONE / 2, pm)
    }
    /// `k` random sequences × `shots` each at length `m` — measured survival.
    pub fn measure(&self, m: u32, k: u32, shots: u32) -> i64 {
        let s = self.survival(m);
        let mut st = self.seed.wrapping_mul(0x27D4_EB2F).wrapping_add(m as u64);
        let mut hits = 0u64;
        let total = (k.max(1) as u64) * shots.max(1) as u64;
        for _ in 0..total {
            let r = (splitmix64(&mut st) >> (64 - CAL_FRAC)) as i64;
            if r < s {
                hits += 1;
            }
        }
        ((hits as i128 * CAL_ONE as i128) / total as i128) as i64
    }
}

#[derive(Clone, Debug, PartialEq, Eq)]
pub struct RbResult {
    /// Fitted depolarizing parameter p.
    pub p_est_fx: i64,
    /// Average 1-qubit gate fidelity `F = 1 − (1−p)/2`.
    pub fidelity_fx: i64,
    pub curve: Vec<(u32, i64)>,
}

impl RbResult {
    pub fn hash(&self) -> [u8; 32] {
        let mut h = blake3::Hasher::new();
        h.update(b"wai:quantum-cal-rb\x01");
        h.update(&self.p_est_fx.to_le_bytes());
        h.update(&self.fidelity_fx.to_le_bytes());
        for (m, s) in &self.curve {
            h.update(&m.to_le_bytes());
            h.update(&s.to_le_bytes());
        }
        *h.finalize().as_bytes()
    }
    pub fn artifacts(&self) -> CalArtifacts {
        let mut ev = Vec::new();
        for (m, s) in &self.curve {
            ev.extend_from_slice(&m.to_le_bytes());
            ev.extend_from_slice(&s.to_le_bytes());
        }
        CalArtifacts {
            target: "1q_gate".into(),
            config_bytes: self.p_est_fx.to_le_bytes().to_vec(),
            evidence_kind: "randomized_benchmarking".into(),
            evidence_bytes: ev,
            summary: format!("p={} fidelity={} (2^{CAL_FRAC})", self.p_est_fx, self.fidelity_fx),
        }
    }
}

#[derive(Clone, Copy, Debug)]
pub struct RbConfig {
    /// Sequence lengths to sample (inclusive powers/steps chosen by the caller).
    pub lengths: [u32; 8],
    pub k_sequences: u32,
    pub shots: u32,
    pub n_grid: u32,
}

pub fn calibrate_rb(dev: &RbDevice, cfg: &RbConfig) -> RbResult {
    let curve: Vec<(u32, i64)> =
        cfg.lengths.iter().map(|&m| (m, dev.measure(m, cfg.k_sequences, cfg.shots))).collect();
    // fit p: minimize Σ |measured − (½ + ½·p^m)| over candidate p ∈ (0,1)
    let ng = cfg.n_grid.max(8);
    let mut best_p = CAL_ONE;
    let mut best_err = i128::MAX;
    for g in 1..ng {
        let p = ((CAL_ONE as i128 * g as i128) / ng as i128) as i64;
        let mut err = 0i128;
        for &(m, s) in &curve {
            let mut pm = CAL_ONE;
            for _ in 0..m {
                pm = fmul(pm, p);
            }
            let pred = CAL_ONE / 2 + fmul(CAL_ONE / 2, pm);
            err += (s - pred).unsigned_abs() as i128;
        }
        if err < best_err {
            best_err = err;
            best_p = p;
        }
    }
    let fidelity = CAL_ONE - (CAL_ONE - best_p) / 2;
    RbResult { p_est_fx: best_p, fidelity_fx: fidelity, curve }
}

// ===========================================================================
// Readout discrimination — GMM classification fidelity
// ===========================================================================

/// A dispersive readout projected onto its discrimination axis: `|0⟩ ~ N(s0,σ)`,
/// `|1⟩ ~ N(s1,σ)`. The calibrator gets labeled shots and finds the threshold.
#[derive(Clone, Copy, Debug)]
pub struct ReadoutDevice {
    pub s0_fx: i64,
    pub s1_fx: i64,
    pub sigma_fx: i64,
    pub seed: u64,
}

impl ReadoutDevice {
    fn gaussian(mean: i64, sigma: i64, st: &mut u64) -> i64 {
        // sum of 12 U[0,1) − 6 ≈ N(0,1)
        let mut acc = 0i64;
        for _ in 0..12 {
            acc += (splitmix64(st) >> (64 - CAL_FRAC)) as i64;
        }
        let z = acc - 6 * CAL_ONE;
        mean + fmul(z, sigma)
    }
    /// `n` labeled shots per class: `(zeros, ones)` readout values.
    pub fn dataset(&self, n: u32) -> (Vec<i64>, Vec<i64>) {
        let mut st0 = self.seed ^ 0xA5A5;
        let mut st1 = self.seed ^ 0x5A5A;
        let zeros = (0..n).map(|_| Self::gaussian(self.s0_fx, self.sigma_fx, &mut st0)).collect();
        let ones = (0..n).map(|_| Self::gaussian(self.s1_fx, self.sigma_fx, &mut st1)).collect();
        (zeros, ones)
    }
}

#[derive(Clone, Debug, PartialEq, Eq)]
pub struct ReadoutResult {
    pub threshold_fx: i64,
    pub mean0_fx: i64,
    pub mean1_fx: i64,
    /// Assignment fidelity `1 − ½(P(1|0)+P(0|1))` at the fitted threshold.
    pub fidelity_fx: i64,
    pub n_shots: u32,
}

impl ReadoutResult {
    pub fn hash(&self) -> [u8; 32] {
        let mut h = blake3::Hasher::new();
        h.update(b"wai:quantum-cal-readout\x01");
        for v in [self.threshold_fx, self.mean0_fx, self.mean1_fx, self.fidelity_fx] {
            h.update(&v.to_le_bytes());
        }
        h.update(&self.n_shots.to_le_bytes());
        *h.finalize().as_bytes()
    }
    pub fn artifacts(&self) -> CalArtifacts {
        let mut cfg = Vec::new();
        cfg.extend_from_slice(&self.threshold_fx.to_le_bytes());
        CalArtifacts {
            target: "readout".into(),
            config_bytes: cfg,
            evidence_kind: "readout_gmm".into(),
            evidence_bytes: {
                let mut e = Vec::new();
                for v in [self.mean0_fx, self.mean1_fx, self.fidelity_fx] {
                    e.extend_from_slice(&v.to_le_bytes());
                }
                e.extend_from_slice(&self.n_shots.to_le_bytes());
                e
            },
            summary: format!("threshold={} fidelity={} (2^{CAL_FRAC})", self.threshold_fx, self.fidelity_fx),
        }
    }
}

pub fn calibrate_readout(dev: &ReadoutDevice, n_per_class: u32) -> ReadoutResult {
    let (zeros, ones) = dev.dataset(n_per_class);
    let mean = |v: &[i64]| (v.iter().map(|&x| x as i128).sum::<i128>() / v.len().max(1) as i128) as i64;
    let m0 = mean(&zeros);
    let m1 = mean(&ones);
    let threshold = (m0 + m1) / 2;
    // classify: predict 1 if value is on the |1⟩ side of the threshold
    let high1 = m1 >= m0;
    let mis0 = zeros.iter().filter(|&&x| (x >= threshold) == high1).count();
    let mis1 = ones.iter().filter(|&&x| (x >= threshold) != high1).count();
    let n = n_per_class.max(1) as i128;
    let err = ((mis0 as i128 * CAL_ONE as i128 / n) + (mis1 as i128 * CAL_ONE as i128 / n)) / 2;
    ReadoutResult {
        threshold_fx: threshold,
        mean0_fx: m0,
        mean1_fx: m1,
        fidelity_fx: CAL_ONE - err as i64,
        n_shots: n_per_class,
    }
}

// ===========================================================================
// Charge-stability tuning — silicon-spin double-dot charge-state autotuning:
// find the gate-voltage
// operating point for a target charge configuration from a 2D charge-stability
// scan. The honeycomb is a clean integer function of the *hidden* capacitance
// model; the tuner discovers the operating point from the scan, byte-exactly —
// exactly what a silicon-spin autotuner does, minus the float and plus a receipt.
// ===========================================================================

/// A double quantum dot under a constant-interaction capacitance model. Two
/// plunger gates set the charge configuration `(n_L, n_R)`; cross-capacitance
/// (`beta`) tilts the transitions into the classic honeycomb. The lever arms and
/// offsets are the hidden truth the tuner never sees.
#[derive(Clone, Copy, Debug)]
pub struct ChargeDevice {
    pub alpha_l_fx: i64,
    pub alpha_r_fx: i64,
    pub beta_lr_fx: i64,
    pub beta_rl_fx: i64,
    pub off_l_fx: i64,
    pub off_r_fx: i64,
}

impl ChargeDevice {
    /// The ground-state charge configuration `(n_L, n_R)` at gate voltages
    /// `(v_l, v_r)` — what the charge sensor reads. Pure integer; the `>> CAL_FRAC`
    /// is a floor (arithmetic shift), clamped to non-negative occupation.
    pub fn charge_state(&self, v_l: i64, v_r: i64) -> (i32, i32) {
        let e_l = fmul(self.alpha_l_fx, v_l) + fmul(self.beta_lr_fx, v_r) + self.off_l_fx;
        let e_r = fmul(self.alpha_r_fx, v_r) + fmul(self.beta_rl_fx, v_l) + self.off_r_fx;
        ((e_l >> CAL_FRAC).max(0) as i32, (e_r >> CAL_FRAC).max(0) as i32)
    }
}

/// Scan window and resolution for [`tune_charge_state`].
#[derive(Clone, Copy, Debug)]
pub struct ChargeScanConfig {
    pub v_max_fx: i64,
    pub n_grid: u32,
}
impl Default for ChargeScanConfig {
    fn default() -> Self {
        ChargeScanConfig { v_max_fx: 8 * CAL_ONE, n_grid: 96 }
    }
}

/// The result of an autonomous charge-state tune-up. Deterministic given
/// `(device, target, cfg)`.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct ChargeResult {
    /// The gate-voltage operating point — the centroid of the target charge cell.
    pub v_l_op_fx: i64,
    pub v_r_op_fx: i64,
    pub target: (i32, i32),
    pub found: bool,
    pub n_grid: u32,
    pub v_max_fx: i64,
    /// Row-major `n_grid × n_grid` total-charge signal `n_L + n_R` (the honeycomb).
    pub signal: Vec<i32>,
}

impl ChargeResult {
    pub fn hash(&self) -> [u8; 32] {
        let mut h = blake3::Hasher::new();
        h.update(b"wai:quantum-cal-chargestab\x01");
        h.update(&self.v_l_op_fx.to_le_bytes());
        h.update(&self.v_r_op_fx.to_le_bytes());
        h.update(&[self.target.0 as u8, self.target.1 as u8, self.found as u8]);
        for s in &self.signal {
            h.update(&s.to_le_bytes());
        }
        *h.finalize().as_bytes()
    }
    pub fn artifacts(&self) -> CalArtifacts {
        let mut cfg = Vec::new();
        cfg.extend_from_slice(&self.v_l_op_fx.to_le_bytes());
        cfg.extend_from_slice(&self.v_r_op_fx.to_le_bytes());
        let mut ev = Vec::with_capacity(self.signal.len() * 4);
        for s in &self.signal {
            ev.extend_from_slice(&s.to_le_bytes());
        }
        CalArtifacts {
            target: format!("charge_state({},{})", self.target.0, self.target.1),
            config_bytes: cfg,
            evidence_kind: "charge_stability_scan".into(),
            evidence_bytes: ev,
            summary: format!(
                "op=({},{}) target=({},{}) found={} (2^{CAL_FRAC})",
                self.v_l_op_fx, self.v_r_op_fx, self.target.0, self.target.1, self.found
            ),
        }
    }
}

/// Autonomously tune the double dot to a target charge configuration: scan the
/// charge-stability diagram over `[0, v_max]²`, identify every grid point in the
/// target cell, and return the operating point — its centroid (which lies inside
/// the convex cell, so re-reading it yields the target). Pure integer,
/// byte-identical everywhere. The tuner is blind to the capacitance model.
pub fn tune_charge_state(
    dev: &ChargeDevice,
    target: (i32, i32),
    cfg: &ChargeScanConfig,
) -> ChargeResult {
    let n = cfg.n_grid.max(2);
    let mut signal = Vec::with_capacity((n * n) as usize);
    let (mut sum_l, mut sum_r, mut count) = (0i128, 0i128, 0i128);
    for iy in 0..n {
        let v_r = ((cfg.v_max_fx as i128 * iy as i128) / (n - 1) as i128) as i64;
        for ix in 0..n {
            let v_l = ((cfg.v_max_fx as i128 * ix as i128) / (n - 1) as i128) as i64;
            let (nl, nr) = dev.charge_state(v_l, v_r);
            signal.push(nl + nr);
            if (nl, nr) == target {
                sum_l += v_l as i128;
                sum_r += v_r as i128;
                count += 1;
            }
        }
    }
    let found = count > 0;
    let (v_l_op_fx, v_r_op_fx) =
        if found { ((sum_l / count) as i64, (sum_r / count) as i64) } else { (0, 0) };
    ChargeResult { v_l_op_fx, v_r_op_fx, target, found, n_grid: n, v_max_fx: cfg.v_max_fx, signal }
}

// ===========================================================================
// Two-qubit cross-resonance gate calibration
// (a standard two-qubit tune-up). The target qubit's
// Rabi rate depends on the CONTROL qubit's state; calibrate both conditional
// π-amplitudes, whose difference is the ZX interaction. Model-fit, byte-exact.
// ===========================================================================

/// A cross-resonance pair: the target's π-amplitude differs by the control state
/// (the ZX shift). Both amplitudes are the hidden truth.
#[derive(Clone, Copy, Debug)]
pub struct CrDevice {
    pub a_pi0_fx: i64,
    pub a_pi1_fx: i64,
    pub seed: u64,
}

impl CrDevice {
    /// Measured target excited fraction at amplitude `a` for a given control
    /// state — a conditional Rabi experiment (deterministic shot noise).
    pub fn measure(&self, a_fx: i64, control: bool, shots: u32) -> i64 {
        let a_pi = if control { self.a_pi1_fx } else { self.a_pi0_fx };
        Device { a_pi_fx: a_pi.max(1), seed: self.seed ^ (control as u64) }.measure(a_fx, shots)
    }
}

#[derive(Clone, Copy, Debug)]
pub struct CrConfig {
    pub a_max_fx: i64,
    pub n_points: u32,
    pub shots: u32,
    pub n_fit: u32,
}
impl Default for CrConfig {
    fn default() -> Self {
        CrConfig { a_max_fx: 2 * CAL_ONE, n_points: 41, shots: 2000, n_fit: 41 }
    }
}

#[derive(Clone, Debug, PartialEq, Eq)]
pub struct CrResult {
    pub a_pi0_est_fx: i64,
    pub a_pi1_est_fx: i64,
    /// ZX interaction proxy: `|a_pi1 − a_pi0|` (the control-conditioned shift).
    pub zx_shift_fx: i64,
    pub sweep0: Vec<(i64, i64)>,
    pub sweep1: Vec<(i64, i64)>,
}

impl CrResult {
    pub fn hash(&self) -> [u8; 32] {
        let mut h = blake3::Hasher::new();
        h.update(b"wai:quantum-cal-cr\x01");
        for v in [self.a_pi0_est_fx, self.a_pi1_est_fx, self.zx_shift_fx] {
            h.update(&v.to_le_bytes());
        }
        for (a, p) in self.sweep0.iter().chain(self.sweep1.iter()) {
            h.update(&a.to_le_bytes());
            h.update(&p.to_le_bytes());
        }
        *h.finalize().as_bytes()
    }
    pub fn artifacts(&self) -> CalArtifacts {
        let mut cfg = Vec::new();
        cfg.extend_from_slice(&self.a_pi0_est_fx.to_le_bytes());
        cfg.extend_from_slice(&self.a_pi1_est_fx.to_le_bytes());
        let mut ev = Vec::new();
        for (a, p) in self.sweep0.iter().chain(self.sweep1.iter()) {
            ev.extend_from_slice(&a.to_le_bytes());
            ev.extend_from_slice(&p.to_le_bytes());
        }
        CalArtifacts {
            target: "cz(control,target)".into(),
            config_bytes: cfg,
            evidence_kind: "cross_resonance_sweep".into(),
            evidence_bytes: ev,
            summary: format!(
                "a_pi0={} a_pi1={} zx_shift={} (2^{CAL_FRAC})",
                self.a_pi0_est_fx, self.a_pi1_est_fx, self.zx_shift_fx
            ),
        }
    }
}

/// Fit a Rabi π-amplitude by whole-curve model matching (shared with the
/// single-qubit routine's approach): minimize `Σ|measured − sin²(π·a/2a_π)|`.
fn fit_pi_amplitude(sweep: &[(i64, i64)], a_max: i64, coarse: u32, fine: u32) -> i64 {
    let fit = |lo: i64, hi: i64, steps: u32| -> i64 {
        let steps = steps.max(2);
        let mut best = lo.max(1);
        let mut best_err = i128::MAX;
        for g in 1..=steps {
            let a_pi = lo + ((hi - lo) as i128 * g as i128 / steps as i128) as i64;
            if a_pi <= 0 {
                continue;
            }
            let mut err = 0i128;
            for &(a, p) in sweep {
                let theta = ((PI_FX as i128 * a as i128) / (2 * a_pi as i128)) as i64;
                err += (p - sin2_fx(theta)).unsigned_abs() as i128;
            }
            if err < best_err {
                best_err = err;
                best = a_pi;
            }
        }
        best
    };
    let a_coarse = fit(a_max / coarse.max(2) as i64, a_max, coarse.max(24));
    let step = a_max / coarse.max(2) as i64;
    fit((a_coarse - step).max(1), a_coarse + step, fine.max(24))
}

/// Calibrate a cross-resonance two-qubit gate: sweep the drive amplitude with the
/// control prepared in `|0⟩` and `|1⟩`, model-fit both conditional π-amplitudes,
/// and report their difference (the ZX interaction). Pure integer, byte-exact.
pub fn calibrate_cr(dev: &CrDevice, cfg: &CrConfig) -> CrResult {
    let n = cfg.n_points.max(4);
    let (mut sweep0, mut sweep1) = (Vec::with_capacity(n as usize), Vec::with_capacity(n as usize));
    for i in 0..n {
        let a = ((cfg.a_max_fx as i128 * i as i128) / (n - 1) as i128) as i64;
        sweep0.push((a, dev.measure(a, false, cfg.shots)));
        sweep1.push((a, dev.measure(a, true, cfg.shots)));
    }
    let a0 = fit_pi_amplitude(&sweep0, cfg.a_max_fx, cfg.n_points, cfg.n_fit);
    let a1 = fit_pi_amplitude(&sweep1, cfg.a_max_fx, cfg.n_points, cfg.n_fit);
    CrResult { a_pi0_est_fx: a0, a_pi1_est_fx: a1, zx_shift_fx: (a1 - a0).abs(), sweep0, sweep1 }
}

// ===========================================================================
// Crosstalk characterization — the N×N drive-response matrix
// (crosstalk mitigation, step one: measure the matrix).
// ===========================================================================

/// An `n`-qubit device with a hidden crosstalk matrix: driving qubit `i` at
/// amplitude `a` rotates qubit `j` by `matrix[i*n+j]·a` (diagonal ≈ 1, off-
/// diagonal small). Row-major, fixed-point.
#[derive(Clone, Debug)]
pub struct CrosstalkDevice {
    pub n: usize,
    pub matrix_fx: Vec<i64>,
    pub seed: u64,
}

impl CrosstalkDevice {
    /// Measured excited fraction of qubit `j` when driving qubit `i` at amplitude
    /// `a`: `P = sin²(matrix[i,j]·a / 2)`.
    pub fn measure(&self, i: usize, j: usize, a_fx: i64, shots: u32) -> i64 {
        let theta = fmul(self.matrix_fx[i * self.n + j], a_fx);
        let p = sin2_fx(theta / 2);
        let mut st = self.seed.wrapping_mul(0x1000_0001).wrapping_add((i * 97 + j) as u64);
        let mut hits = 0u64;
        for _ in 0..shots {
            if ((splitmix64(&mut st) >> (64 - CAL_FRAC)) as i64) < p {
                hits += 1;
            }
        }
        ((hits as i128 * CAL_ONE as i128) / shots.max(1) as i128) as i64
    }
}

#[derive(Clone, Debug, PartialEq, Eq)]
pub struct CrosstalkResult {
    pub n: usize,
    /// Recovered `n×n` crosstalk coefficients (row-major, fixed-point).
    pub matrix_fx: Vec<i64>,
    /// Worst off-diagonal coefficient (the crosstalk a mitigation must cancel).
    pub worst_offdiag_fx: i64,
}

impl CrosstalkResult {
    pub fn hash(&self) -> [u8; 32] {
        let mut h = blake3::Hasher::new();
        h.update(b"wai:quantum-cal-crosstalk\x01");
        h.update(&(self.n as u64).to_le_bytes());
        for v in &self.matrix_fx {
            h.update(&v.to_le_bytes());
        }
        *h.finalize().as_bytes()
    }
    pub fn artifacts(&self) -> CalArtifacts {
        let mut ev = Vec::new();
        for v in &self.matrix_fx {
            ev.extend_from_slice(&v.to_le_bytes());
        }
        CalArtifacts {
            target: format!("crosstalk_{}x{}", self.n, self.n),
            config_bytes: self.worst_offdiag_fx.to_le_bytes().to_vec(),
            evidence_kind: "crosstalk_matrix".into(),
            evidence_bytes: ev,
            summary: format!("n={} worst_offdiag={} (2^{CAL_FRAC})", self.n, self.worst_offdiag_fx),
        }
    }
}

/// Characterize the crosstalk matrix: for each `(drive i, measure j)`, sweep a
/// small amplitude, and recover `matrix[i,j]` from the initial response slope
/// (`P ≈ (matrix·a/2)²` for small angle ⇒ `matrix ≈ 2·√P / a`). Pure integer.
pub fn characterize_crosstalk(dev: &CrosstalkDevice, probe_a_fx: i64, shots: u32) -> CrosstalkResult {
    let n = dev.n;
    let mut matrix = vec![0i64; n * n];
    // fit each coefficient by the whole-response model over a small sweep
    let sweeps = 12u32;
    for i in 0..n {
        for j in 0..n {
            let mut sweep = Vec::with_capacity(sweeps as usize);
            for s in 1..=sweeps {
                let a = (probe_a_fx as i128 * s as i128 / sweeps as i128) as i64;
                sweep.push((a, dev.measure(i, j, a, shots)));
            }
            // search coefficient c minimizing Σ|measured − sin²(c·a/2)|
            let mut best = 0i64;
            let mut best_err = i128::MAX;
            let c_max = 2 * CAL_ONE;
            for g in 0..=200 {
                let c = (c_max as i128 * g as i128 / 200) as i64;
                let mut err = 0i128;
                for &(a, p) in &sweep {
                    err += (p - sin2_fx(fmul(c, a) / 2)).unsigned_abs() as i128;
                }
                if err < best_err {
                    best_err = err;
                    best = c;
                }
            }
            matrix[i * n + j] = best;
        }
    }
    let worst = (0..n)
        .flat_map(|i| (0..n).filter(move |&j| i != j).map(move |j| (i, j)))
        .map(|(i, j)| matrix[i * n + j].abs())
        .max()
        .unwrap_or(0);
    CrosstalkResult { n, matrix_fx: matrix, worst_offdiag_fx: worst }
}

// ===========================================================================
// DRAG / leakage pulse calibration — the pulse-shaping (optimal-control) primitive
// (leakage suppression). Optimize the DRAG coefficient β that
// minimizes leakage to |2⟩. A parabolic minimization — the simplest honest
// representative of the control-optimization layer (NOT full robust control).
// ===========================================================================

/// A qubit whose leakage-to-`|2⟩` is minimized at a hidden DRAG coefficient
/// `beta_opt` (a parabola in β). `curvature` sets how sharply leakage rises.
#[derive(Clone, Copy, Debug)]
pub struct DragDevice {
    pub beta_opt_fx: i64,
    pub curvature_fx: i64,
    pub seed: u64,
}

impl DragDevice {
    /// Measured leakage signal at DRAG coefficient `β`: `L = curvature·(β−β_opt)²`
    /// (clamped to `[0,1]`), with deterministic measurement noise.
    pub fn measure_leakage(&self, beta_fx: i64, shots: u32) -> i64 {
        let d = beta_fx - self.beta_opt_fx;
        let l = fmul(self.curvature_fx, fmul(d, d)).clamp(0, CAL_ONE);
        let mut st = self.seed.wrapping_mul(0x9E37).wrapping_add(beta_fx as u64);
        let mut hits = 0u64;
        for _ in 0..shots {
            if ((splitmix64(&mut st) >> (64 - CAL_FRAC)) as i64) < l {
                hits += 1;
            }
        }
        ((hits as i128 * CAL_ONE as i128) / shots.max(1) as i128) as i64
    }
}

#[derive(Clone, Debug, PartialEq, Eq)]
pub struct DragResult {
    pub beta_opt_est_fx: i64,
    /// Residual leakage at the optimized β (`≈ 0` for a well-shaped pulse).
    pub residual_leakage_fx: i64,
    pub sweep: Vec<(i64, i64)>,
}

impl DragResult {
    pub fn hash(&self) -> [u8; 32] {
        let mut h = blake3::Hasher::new();
        h.update(b"wai:quantum-cal-drag\x01");
        h.update(&self.beta_opt_est_fx.to_le_bytes());
        h.update(&self.residual_leakage_fx.to_le_bytes());
        for (b, l) in &self.sweep {
            h.update(&b.to_le_bytes());
            h.update(&l.to_le_bytes());
        }
        *h.finalize().as_bytes()
    }
    pub fn artifacts(&self) -> CalArtifacts {
        let mut ev = Vec::new();
        for (b, l) in &self.sweep {
            ev.extend_from_slice(&b.to_le_bytes());
            ev.extend_from_slice(&l.to_le_bytes());
        }
        CalArtifacts {
            target: "drag_beta".into(),
            config_bytes: self.beta_opt_est_fx.to_le_bytes().to_vec(),
            evidence_kind: "drag_leakage_sweep".into(),
            evidence_bytes: ev,
            summary: format!(
                "beta_opt={} residual_leakage={} (2^{CAL_FRAC})",
                self.beta_opt_est_fx, self.residual_leakage_fx
            ),
        }
    }
}

/// Optimize the DRAG coefficient: sweep β over `[beta_min, beta_max]`, measure
/// leakage, and take the minimum (parabolic minimization — the pulse-shaping
/// optimization that suppresses leakage). Pure integer, byte-exact.
pub fn calibrate_drag(dev: &DragDevice, beta_min_fx: i64, beta_max_fx: i64, n: u32, shots: u32) -> DragResult {
    let n = n.max(4);
    let mut sweep = Vec::with_capacity(n as usize);
    let (mut best_b, mut best_l) = (beta_min_fx, i64::MAX);
    for i in 0..n {
        let b = beta_min_fx + ((beta_max_fx - beta_min_fx) as i128 * i as i128 / (n - 1) as i128) as i64;
        let l = dev.measure_leakage(b, shots);
        sweep.push((b, l));
        if l < best_l {
            best_l = l;
            best_b = b;
        }
    }
    // parabolic refinement around the discrete minimum
    let step = (beta_max_fx - beta_min_fx) / (n as i64 - 1);
    let refine = |lo: i64, hi: i64| -> (i64, i64) {
        let (mut bb, mut bl) = (lo, i64::MAX);
        for g in 0..=60 {
            let b = lo + ((hi - lo) as i128 * g as i128 / 60) as i64;
            let l = dev.measure_leakage(b, shots);
            if l < bl {
                bl = l;
                bb = b;
            }
        }
        (bb, bl)
    };
    let (beta, res) = refine((best_b - step).max(beta_min_fx), (best_b + step).min(beta_max_fx));
    DragResult { beta_opt_est_fx: beta, residual_leakage_fx: res, sweep }
}

// ===========================================================================
// Robust optimal control — composite-pulse design.
// Numerically OPTIMIZE a control pulse to implement a π-gate that is robust to
// amplitude noise: minimize the worst-case gate error over an amplitude-error
// range. Bloch-sphere evolution in fixed-point, no float — the optimizer
// discovers the noise-robust composite (vs a fragile naive square pulse).
// ===========================================================================

/// Rotate a Bloch vector `v` about the XY-plane axis `(cos φ, sin φ, 0)` by angle
/// `alpha` (Rodrigues' formula), fixed-point. The rotation axis phase `phi` is
/// the pulse phase; `alpha` the rotation angle (∝ amplitude·time).
fn bloch_rotate(v: [i64; 3], phi: i64, alpha: i64) -> [i64; 3] {
    let (n0, n1) = (cos_fx(phi), sin_fx(phi));
    let (c, s) = (cos_fx(alpha), sin_fx(alpha));
    let ndotv = fmul(n0, v[0]) + fmul(n1, v[1]); // n2 = 0
    let cross = [fmul(n1, v[2]), -fmul(n0, v[2]), fmul(n0, v[1]) - fmul(n1, v[0])];
    let n = [n0, n1, 0];
    let omc = CAL_ONE - c;
    let mut out = [0i64; 3];
    for k in 0..3 {
        out[k] = fmul(v[k], c) + fmul(cross[k], s) + fmul(n[k], fmul(ndotv, omc));
    }
    out
}

/// Apply a composite pulse `[(phase, angle)…]` to `|0⟩` (Bloch `+z`) under a
/// fractional amplitude error `eps` (every angle scaled by `1+eps`), returning
/// the excited-state population `(1 − z)/2` — the population-transfer fidelity of
/// a π-gate. Clamped to `[0, 1]`.
pub fn apply_composite(segments: &[(i64, i64)], eps_fx: i64) -> i64 {
    let mut v = [0i64, 0, CAL_ONE];
    for &(phi, alpha) in segments {
        let a = alpha + fmul(alpha, eps_fx);
        v = bloch_rotate(v, phi, a);
    }
    ((CAL_ONE - v[2]) / 2).clamp(0, CAL_ONE)
}

#[derive(Clone, Debug, PartialEq, Eq)]
pub struct RobustResult {
    /// Optimized symmetric composite `[(0,a1), (phi,a2), (0,a1)]`.
    pub a1_fx: i64,
    pub a2_fx: i64,
    pub phi_fx: i64,
    pub eps_axis_fx: Vec<i64>,
    /// Gate error `1 − P` vs `eps` for the naive single π-pulse.
    pub naive_infid_fx: Vec<i64>,
    /// Gate error vs `eps` for the optimized robust pulse.
    pub robust_infid_fx: Vec<i64>,
    pub worst_naive_fx: i64,
    pub worst_robust_fx: i64,
}

impl RobustResult {
    pub fn hash(&self) -> [u8; 32] {
        let mut h = blake3::Hasher::new();
        h.update(b"wai:quantum-cal-robust\x01");
        for v in [self.a1_fx, self.a2_fx, self.phi_fx, self.worst_naive_fx, self.worst_robust_fx] {
            h.update(&v.to_le_bytes());
        }
        for v in self.naive_infid_fx.iter().chain(self.robust_infid_fx.iter()) {
            h.update(&v.to_le_bytes());
        }
        *h.finalize().as_bytes()
    }
    pub fn artifacts(&self) -> CalArtifacts {
        let mut cfg = Vec::new();
        for v in [self.a1_fx, self.a2_fx, self.phi_fx] {
            cfg.extend_from_slice(&v.to_le_bytes());
        }
        let mut ev = Vec::new();
        for v in self.naive_infid_fx.iter().chain(self.robust_infid_fx.iter()) {
            ev.extend_from_slice(&v.to_le_bytes());
        }
        CalArtifacts {
            target: "x_pi_robust".into(),
            config_bytes: cfg,
            evidence_kind: "robust_pulse_infidelity".into(),
            evidence_bytes: ev,
            summary: format!(
                "worst_naive={} worst_robust={} (2^{CAL_FRAC})",
                self.worst_naive_fx, self.worst_robust_fx
            ),
        }
    }
}

/// Optimize a robust π-gate: search a symmetric 3-segment composite
/// `[(0,a1),(phi,a2),(0,a1)]` for the one that is a valid π-gate at nominal
/// amplitude AND minimizes the worst-case gate error over `eps ∈ [−err_max,
/// err_max]`. Pure integer — the optimizer's choice is byte-identical everywhere.
pub fn optimize_robust_pulse(err_max_fx: i64, n_eps: u32) -> RobustResult {
    let n_eps = n_eps.max(3);
    let eps_axis: Vec<i64> = (0..n_eps)
        .map(|i| -err_max_fx + (2 * err_max_fx as i128 * i as i128 / (n_eps - 1) as i128) as i64)
        .collect();
    let worst = |segs: &[(i64, i64)]| -> (Vec<i64>, i64) {
        let mut curve = Vec::with_capacity(eps_axis.len());
        let mut w = 0i64;
        for &e in &eps_axis {
            let inf = CAL_ONE - apply_composite(segs, e);
            w = w.max(inf);
            curve.push(inf);
        }
        (curve, w)
    };

    let (naive_curve, worst_naive) = worst(&[(0, PI_FX)]);

    let (mut best, mut best_w) = ((PI_FX / 2, PI_FX, 0i64), i64::MAX);
    let a1_lo = PI_FX / 4;
    let quarter = PI_FX / 2;
    for i in 0..=20 {
        let a1 = a1_lo + ((PI_FX - a1_lo) as i128 * i as i128 / 20) as i64;
        for j in 0..=20 {
            let a2 = quarter + ((3 * PI_FX / 2 - quarter) as i128 * j as i128 / 20) as i64;
            for k in 0..=24 {
                let phi = (TWO_PI_FX as i128 * k as i128 / 24) as i64;
                let segs = [(0, a1), (phi, a2), (0, a1)];
                if (CAL_ONE - apply_composite(&segs, 0)).abs() > CAL_ONE / 100 {
                    continue; // must be a valid π-gate at nominal amplitude
                }
                let (_, w) = worst(&segs);
                if w < best_w {
                    best_w = w;
                    best = (a1, a2, phi);
                }
            }
        }
    }
    let best_segs = [(0, best.0), (best.2, best.1), (0, best.0)];
    let (robust_curve, worst_robust) = worst(&best_segs);
    RobustResult {
        a1_fx: best.0,
        a2_fx: best.1,
        phi_fx: best.2,
        eps_axis_fx: eps_axis,
        naive_infid_fx: naive_curve,
        robust_infid_fx: robust_curve,
        worst_naive_fx: worst_naive,
        worst_robust_fx: worst_robust,
    }
}

// ===========================================================================
// Generic seal — any routine's artifacts → wai.quantum.calibration receipt
// ===========================================================================

#[cfg(feature = "quantum_ops")]
mod seal_any {
    use super::CalArtifacts;
    use crate::quantum_ops::{content_hash, CalibrationReceipt, Evidence, GrantRef};
    use ed25519_dalek::SigningKey;

    /// Seal any calibration routine's [`CalArtifacts`] into a receipt: the config
    /// hash pins the calibrated parameter, the evidence content-addresses the
    /// dataset, the joules are the run's measured energy, the grant authorizes it.
    #[allow(clippy::too_many_arguments)]
    pub fn seal_artifacts(
        signer: &SigningKey,
        signer_id: impl Into<String>,
        device_id: impl Into<String>,
        art: &CalArtifacts,
        joules_micro: u64,
        grant: GrantRef,
        parent: Option<[u8; 32]>,
    ) -> CalibrationReceipt {
        let ev = vec![Evidence {
            kind: art.evidence_kind.clone(),
            blob_hash: content_hash(&art.evidence_bytes),
            summary: art.summary.clone(),
        }];
        CalibrationReceipt::seal(
            signer, signer_id, device_id, art.target.clone(),
            content_hash(&art.config_bytes), ev, joules_micro, grant, parent,
        )
    }
}

#[cfg(feature = "quantum_ops")]
pub use seal_any::seal_artifacts;

// ===========================================================================
// Calibration-graph orchestration — a bring-up campaign as a signed lineage
// (a node graph of calibrations; here each node's receipt CHAINS to the previous
// via parent_receipt_hash, so the whole campaign is one portable, verifiable lineage).
// ===========================================================================

#[cfg(feature = "quantum_ops")]
mod campaign {
    use super::*;
    use crate::quantum_ops::{CalibrationReceipt, GrantRef};
    use ed25519_dalek::SigningKey;

    /// One node of a calibration campaign — a routine plus its inputs. `run`
    /// dispatches to the deterministic engine and returns its receipt artifacts.
    #[allow(clippy::large_enum_variant)]
    #[derive(Clone, Debug)]
    pub enum CalStep {
        Rabi(Device, TuneConfig),
        Ramsey(RamseyDevice, RamseyConfig),
        Rb(RbDevice, RbConfig),
        Readout(ReadoutDevice, u32),
        Charge(ChargeDevice, (i32, i32), ChargeScanConfig),
    }

    impl CalStep {
        fn run(&self) -> (CalArtifacts, &'static str) {
            match self {
                CalStep::Rabi(d, c) => (calibrate_rabi(d, c).artifacts(), "rabi"),
                CalStep::Ramsey(d, c) => (calibrate_ramsey(d, c).artifacts(), "ramsey"),
                CalStep::Rb(d, c) => (calibrate_rb(d, c).artifacts(), "rb"),
                CalStep::Readout(d, n) => (calibrate_readout(d, *n).artifacts(), "readout"),
                CalStep::Charge(d, t, c) => (tune_charge_state(d, *t, c).artifacts(), "charge"),
            }
        }
    }

    /// A completed calibration campaign: the ordered receipts (each chained to the
    /// one before) and their routine kinds.
    #[derive(Clone, Debug)]
    pub struct Campaign {
        pub receipts: Vec<CalibrationReceipt>,
        pub kinds: Vec<String>,
    }

    impl Campaign {
        /// A single portable id for the whole campaign — BLAKE3 over the ordered
        /// per-node receipt hashes.
        pub fn campaign_hash(&self) -> [u8; 32] {
            let mut h = blake3::Hasher::new();
            h.update(b"wai:quantum-cal-campaign\x01");
            for r in &self.receipts {
                h.update(&r.receipt_hash());
            }
            *h.finalize().as_bytes()
        }

        /// The whole lineage is intact: every receipt verifies, and each node's
        /// `parent_receipt_hash` is exactly the previous node's `receipt_hash`
        /// (the first is a root). A third party confirms the bring-up ran in this
        /// order, each step authorized + metered — a signed calibration campaign.
        pub fn verify_chain(&self) -> bool {
            let mut prev: Option<[u8; 32]> = None;
            for r in &self.receipts {
                if !r.verify() || r.parent_receipt_hash != prev {
                    return false;
                }
                prev = Some(r.receipt_hash());
            }
            true
        }
    }

    /// Run a calibration campaign: execute each step in order, sealing a
    /// `wai.quantum.calibration` receipt chained to the previous step. Every
    /// engine is deterministic, so the whole campaign — receipts and
    /// `campaign_hash` — is byte-identical on every machine.
    pub fn run_campaign(
        signer: &SigningKey,
        signer_id: &str,
        device_id: &str,
        joules_each: u64,
        grant: &GrantRef,
        steps: &[CalStep],
    ) -> Campaign {
        let mut receipts = Vec::with_capacity(steps.len());
        let mut kinds = Vec::with_capacity(steps.len());
        let mut parent: Option<[u8; 32]> = None;
        for step in steps {
            let (art, kind) = step.run();
            let r = seal_artifacts(signer, signer_id, device_id, &art, joules_each, grant.clone(), parent);
            parent = Some(r.receipt_hash());
            kinds.push(kind.to_string());
            receipts.push(r);
        }
        Campaign { receipts, kinds }
    }
}

#[cfg(feature = "quantum_ops")]
pub use campaign::{run_campaign, CalStep, Campaign};

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn sine_is_accurate_enough() {
        // sin²(π/2) = 1, sin²(0) = 0, sin²(π/4) = 0.5
        assert_eq!(sin2_fx(0), 0);
        let at_half_pi = sin2_fx(PI_FX / 2);
        assert!((at_half_pi - CAL_ONE).abs() < CAL_ONE / 200, "sin²(π/2)≈1, got {at_half_pi}");
        let at_quarter = sin2_fx(PI_FX / 4);
        assert!((at_quarter - CAL_ONE / 2).abs() < CAL_ONE / 100, "sin²(π/4)≈0.5, got {at_quarter}");
    }

    #[test]
    fn tune_up_finds_the_true_pi_amplitude() {
        // Hidden truth: π-pulse at amplitude 0.80. The calibrator must find it.
        let a_pi_true = (0.80 * CAL_ONE as f64) as i64;
        let device = Device { a_pi_fx: a_pi_true, seed: 0xCAFE };
        let cfg = TuneConfig { a_max_fx: (1.6 * CAL_ONE as f64) as i64, ..Default::default() };
        let r = calibrate_rabi(&device, &cfg);
        // within one fine step (~0.8%) of the hidden truth
        let tol = cfg.a_max_fx / (cfg.n_coarse as i64 - 1) / (cfg.n_fine as i64 - 1) + CAL_ONE / 50;
        assert!(
            (r.a_pi_est_fx - a_pi_true).abs() <= tol,
            "estimate {} vs true {} (tol {tol})",
            r.a_pi_est_fx, a_pi_true
        );
        // a clean π rotation: measured excited population ≈ 1
        assert!(r.peak_prob_fx > CAL_ONE - CAL_ONE / 20, "peak prob {} not ~1", r.peak_prob_fx);
    }

    #[test]
    fn calibration_is_byte_exact_deterministic() {
        let device = Device { a_pi_fx: (0.73 * CAL_ONE as f64) as i64, seed: 7 };
        let cfg = TuneConfig::default();
        let a = calibrate_rabi(&device, &cfg);
        let b = calibrate_rabi(&device, &cfg);
        assert_eq!(a, b);
        assert_eq!(a.hash(), b.hash(), "same (device,cfg) must give the same tune-up hash");
    }

    #[test]
    fn different_devices_tune_differently() {
        let d1 = Device { a_pi_fx: (0.5 * CAL_ONE as f64) as i64, seed: 1 };
        let d2 = Device { a_pi_fx: (1.1 * CAL_ONE as f64) as i64, seed: 1 };
        let cfg = TuneConfig { a_max_fx: 2 * CAL_ONE, ..Default::default() };
        let r1 = calibrate_rabi(&d1, &cfg);
        let r2 = calibrate_rabi(&d2, &cfg);
        assert!(r1.a_pi_est_fx < r2.a_pi_est_fx);
    }

    #[cfg(feature = "quantum_ops")]
    #[test]
    fn calibrate_and_seal_produces_a_verifying_receipt() {
        use crate::quantum_ops::GrantRef;
        use ed25519_dalek::SigningKey;
        let device = Device { a_pi_fx: (0.9 * CAL_ONE as f64) as i64, seed: 0xABCD };
        let cfg = TuneConfig::default();
        let grant = GrantRef {
            grant_hash: [0u8; 32],
            capability: "quantum.calibrate".into(),
            joule_ceiling_micro: 10_000_000,
            funds_ceiling: None,
        };
        let (r, receipt) = calibrate_and_seal(
            &SigningKey::from_bytes(&[9u8; 32]), "did:key:lab", "sim:transmon:q0", "q0",
            &device, &cfg, 3_500_000, grant, None,
        );
        assert!(receipt.verify());
        // the evidence content-addresses the exact sweep the calibrator produced
        assert!(receipt.evidence_matches(0, &sweep_bytes(&r)));
        // and a re-run reproduces that sweep byte-for-byte (the reproducibility claim)
        let r2 = calibrate_rabi(&device, &cfg);
        assert!(receipt.evidence_matches(0, &sweep_bytes(&r2)));
    }

    fn fx(x: f64) -> i64 {
        (x * CAL_ONE as f64) as i64
    }

    #[test]
    fn ramsey_recovers_detuning_and_t2_deterministically() {
        let dev = RamseyDevice { omega_fx: fx(3.0), t2_fx: fx(3.0), seed: 0x1234 };
        let cfg = RamseyConfig {
            tau_max_fx: fx(6.0), n_points: 80, shots: 4000,
            omega_max_fx: fx(6.0), t2_max_fx: fx(6.0), n_grid: 60,
        };
        let r = calibrate_ramsey(&dev, &cfg);
        assert!((r.omega_est_fx - fx(3.0)).abs() < fx(0.4), "omega {} vs 3.0", r.omega_est_fx);
        assert!((r.t2_est_fx - fx(3.0)).abs() < fx(0.8), "t2 {} vs 3.0", r.t2_est_fx);
        // byte-exact determinism
        assert_eq!(r, calibrate_ramsey(&dev, &cfg));
        assert_eq!(r.hash(), calibrate_ramsey(&dev, &cfg).hash());
    }

    #[test]
    fn rb_recovers_gate_fidelity_deterministically() {
        // 1% per-gate error → p = 0.98, F = 0.99.
        let dev = RbDevice { epsilon_fx: fx(0.01), seed: 0xBEEF };
        let cfg = RbConfig {
            lengths: [1, 2, 4, 8, 16, 32, 64, 128], k_sequences: 20, shots: 100, n_grid: 300,
        };
        let r = calibrate_rb(&dev, &cfg);
        assert!((r.p_est_fx - fx(0.98)).abs() < fx(0.02), "p {} vs 0.98", r.p_est_fx);
        assert!((r.fidelity_fx - fx(0.99)).abs() < fx(0.01), "F {} vs 0.99", r.fidelity_fx);
        assert_eq!(r, calibrate_rb(&dev, &cfg));
    }

    #[test]
    fn readout_discriminates_well_separated_blobs_deterministically() {
        // 4σ separation → high assignment fidelity.
        let dev = ReadoutDevice { s0_fx: fx(-2.0), s1_fx: fx(2.0), sigma_fx: fx(1.0), seed: 42 };
        let r = calibrate_readout(&dev, 500);
        assert!(r.fidelity_fx > fx(0.9), "readout fidelity {} not > 0.9", r.fidelity_fx);
        assert!(r.mean0_fx < r.threshold_fx && r.threshold_fx < r.mean1_fx);
        assert_eq!(r, calibrate_readout(&dev, 500));
        assert_eq!(r.hash(), calibrate_readout(&dev, 500).hash());
    }

    #[cfg(feature = "quantum_ops")]
    #[test]
    fn every_routine_seals_a_verifying_receipt() {
        use crate::quantum_ops::GrantRef;
        use ed25519_dalek::SigningKey;
        let signer = SigningKey::from_bytes(&[5u8; 32]);
        let grant = GrantRef::unbounded("quantum.calibrate");

        let ram = calibrate_ramsey(
            &RamseyDevice { omega_fx: fx(2.0), t2_fx: fx(4.0), seed: 1 },
            &RamseyConfig { tau_max_fx: fx(8.0), n_points: 60, shots: 3000, omega_max_fx: fx(6.0), t2_max_fx: fx(8.0), n_grid: 50 },
        );
        let rb = calibrate_rb(
            &RbDevice { epsilon_fx: fx(0.005), seed: 2 },
            &RbConfig { lengths: [1, 2, 4, 8, 16, 32, 64, 128], k_sequences: 10, shots: 100, n_grid: 200 },
        );
        let ro = calibrate_readout(&ReadoutDevice { s0_fx: fx(-1.5), s1_fx: fx(1.5), sigma_fx: fx(1.0), seed: 3 }, 400);

        for art in [ram.artifacts(), rb.artifacts(), ro.artifacts()] {
            let r = seal_artifacts(&signer, "did:key:lab", "sim:transmon:q0", &art, 1_000_000, grant.clone(), None);
            assert!(r.verify(), "receipt for {} must verify", art.evidence_kind);
            assert!(r.evidence_matches(0, &art.evidence_bytes));
        }
    }

    fn dot() -> ChargeDevice {
        ChargeDevice {
            alpha_l_fx: fx(0.4), alpha_r_fx: fx(0.4),
            beta_lr_fx: fx(0.1), beta_rl_fx: fx(0.1),
            off_l_fx: 0, off_r_fx: 0,
        }
    }

    #[test]
    fn charge_tune_is_byte_exact_deterministic() {
        let (d, cfg) = (dot(), ChargeScanConfig::default());
        let a = tune_charge_state(&d, (1, 1), &cfg);
        let b = tune_charge_state(&d, (1, 1), &cfg);
        assert_eq!(a, b);
        assert_eq!(a.hash(), b.hash());
    }

    #[test]
    fn charge_tune_operating_point_reads_the_target() {
        let (d, cfg) = (dot(), ChargeScanConfig::default());
        for target in [(1, 1), (2, 1), (1, 2), (2, 2)] {
            let r = tune_charge_state(&d, target, &cfg);
            assert!(r.found, "target {target:?} not in the scan window");
            // the centroid lies inside the convex cell → re-reading yields the target
            assert_eq!(
                d.charge_state(r.v_l_op_fx, r.v_r_op_fx),
                target,
                "operating point for {target:?} reads the wrong charge state"
            );
        }
    }

    #[test]
    fn charge_tune_higher_target_needs_higher_voltage() {
        let (d, cfg) = (dot(), ChargeScanConfig::default());
        let r1 = tune_charge_state(&d, (1, 1), &cfg);
        let r2 = tune_charge_state(&d, (2, 2), &cfg);
        assert!(r2.v_l_op_fx > r1.v_l_op_fx && r2.v_r_op_fx > r1.v_r_op_fx);
    }

    #[cfg(feature = "quantum_ops")]
    #[test]
    fn charge_tune_seals_a_verifying_receipt() {
        use crate::quantum_ops::GrantRef;
        use ed25519_dalek::SigningKey;
        let r = tune_charge_state(&dot(), (1, 1), &ChargeScanConfig::default());
        let art = r.artifacts();
        let rec = seal_artifacts(
            &SigningKey::from_bytes(&[6u8; 32]), "did:key:lab", "sim:double-dot:d0",
            &art, 2_000_000, GrantRef::unbounded("quantum.calibrate"), None,
        );
        assert!(rec.verify());
        assert!(rec.evidence_matches(0, &art.evidence_bytes));
    }

    #[cfg(feature = "quantum_ops")]
    #[test]
    fn campaign_chains_receipts_into_a_signed_lineage() {
        use crate::quantum_ops::GrantRef;
        use ed25519_dalek::SigningKey;
        let signer = SigningKey::from_bytes(&[7u8; 32]);
        let grant = GrantRef::unbounded("quantum.calibrate");
        let steps = vec![
            CalStep::Readout(
                ReadoutDevice { s0_fx: fx(-1.5), s1_fx: fx(1.5), sigma_fx: fx(1.0), seed: 1 },
                400,
            ),
            CalStep::Rabi(Device { a_pi_fx: fx(0.8), seed: 2 }, TuneConfig::default()),
            CalStep::Charge(dot(), (1, 1), ChargeScanConfig::default()),
        ];
        let camp = run_campaign(&signer, "did:key:lab", "sim:q0", 1_000_000, &grant, &steps);
        assert_eq!(camp.receipts.len(), 3);
        assert_eq!(camp.kinds, ["readout", "rabi", "charge"]);
        assert!(camp.verify_chain(), "the campaign lineage must verify");
        // deterministic: same steps → byte-identical campaign
        let camp2 = run_campaign(&signer, "did:key:lab", "sim:q0", 1_000_000, &grant, &steps);
        assert_eq!(camp.campaign_hash(), camp2.campaign_hash());
    }

    #[cfg(feature = "quantum_ops")]
    #[test]
    fn campaign_tamper_breaks_the_chain() {
        use crate::quantum_ops::GrantRef;
        use ed25519_dalek::SigningKey;
        let steps = vec![
            CalStep::Rabi(Device { a_pi_fx: fx(0.8), seed: 2 }, TuneConfig::default()),
            CalStep::Rb(
                RbDevice { epsilon_fx: fx(0.005), seed: 3 },
                RbConfig { lengths: [1, 2, 4, 8, 16, 32, 64, 128], k_sequences: 8, shots: 80, n_grid: 150 },
            ),
        ];
        let mut camp = run_campaign(
            &SigningKey::from_bytes(&[8u8; 32]), "m", "d", 10,
            &GrantRef::unbounded("quantum.calibrate"), &steps,
        );
        // snap the lineage: re-point the 2nd node's parent
        camp.receipts[1].parent_receipt_hash = Some([0u8; 32]);
        assert!(!camp.verify_chain());
    }

    #[test]
    fn cr_recovers_both_conditional_amplitudes() {
        let dev = CrDevice { a_pi0_fx: fx(0.70), a_pi1_fx: fx(0.95), seed: 11 };
        let cfg = CrConfig { a_max_fx: fx(1.8), ..Default::default() };
        let r = calibrate_cr(&dev, &cfg);
        assert!((r.a_pi0_est_fx - fx(0.70)).abs() < CAL_ONE / 20, "a0 {}", r.a_pi0_est_fx);
        assert!((r.a_pi1_est_fx - fx(0.95)).abs() < CAL_ONE / 20, "a1 {}", r.a_pi1_est_fx);
        assert!(r.zx_shift_fx > 0);
        assert_eq!(r, calibrate_cr(&dev, &cfg));
    }

    #[test]
    fn crosstalk_recovers_the_matrix() {
        let dev = CrosstalkDevice {
            n: 2,
            matrix_fx: vec![fx(1.0), fx(0.12), fx(0.15), fx(1.0)],
            seed: 5,
        };
        let r = characterize_crosstalk(&dev, fx(3.0), 4000);
        for k in 0..4 {
            assert!(
                (r.matrix_fx[k] - dev.matrix_fx[k]).abs() < CAL_ONE / 10,
                "coeff {k}: got {} want {}",
                r.matrix_fx[k], dev.matrix_fx[k]
            );
        }
        assert!(r.worst_offdiag_fx > 0);
        assert_eq!(r, characterize_crosstalk(&dev, fx(3.0), 4000));
    }

    #[test]
    fn drag_finds_the_leakage_minimum() {
        let dev = DragDevice { beta_opt_fx: fx(0.5), curvature_fx: fx(2.0), seed: 9 };
        let r = calibrate_drag(&dev, fx(-1.0), fx(2.0), 41, 3000);
        assert!((r.beta_opt_est_fx - fx(0.5)).abs() < CAL_ONE / 12, "beta {}", r.beta_opt_est_fx);
        assert!(r.residual_leakage_fx < CAL_ONE / 20, "residual {}", r.residual_leakage_fx);
        assert_eq!(r, calibrate_drag(&dev, fx(-1.0), fx(2.0), 41, 3000));
    }

    #[cfg(feature = "quantum_ops")]
    #[test]
    fn new_routines_seal_verifying_receipts() {
        use crate::quantum_ops::GrantRef;
        use ed25519_dalek::SigningKey;
        let signer = SigningKey::from_bytes(&[4u8; 32]);
        let grant = GrantRef::unbounded("quantum.calibrate");
        let cr = calibrate_cr(&CrDevice { a_pi0_fx: fx(0.7), a_pi1_fx: fx(0.95), seed: 1 }, &CrConfig::default());
        let xt = characterize_crosstalk(&CrosstalkDevice { n: 2, matrix_fx: vec![fx(1.0), fx(0.12), fx(0.15), fx(1.0)], seed: 2 }, fx(3.0), 2000);
        let dg = calibrate_drag(&DragDevice { beta_opt_fx: fx(0.5), curvature_fx: fx(2.0), seed: 3 }, fx(-1.0), fx(2.0), 31, 2000);
        for art in [cr.artifacts(), xt.artifacts(), dg.artifacts()] {
            let r = seal_artifacts(&signer, "did:key:lab", "sim:q0", &art, 1_000_000, grant.clone(), None);
            assert!(r.verify(), "receipt for {} must verify", art.evidence_kind);
            assert!(r.evidence_matches(0, &art.evidence_bytes));
        }
    }

    #[test]
    fn sin_cos_are_accurate() {
        assert!(sin_fx(0).abs() < CAL_ONE / 400);
        assert!((sin_fx(PI_FX / 2) - CAL_ONE).abs() < CAL_ONE / 200);
        assert!(sin_fx(PI_FX).abs() < CAL_ONE / 200);
        assert!((cos_fx(0) - CAL_ONE).abs() < CAL_ONE / 200);
        assert!((sin_fx(PI_FX + PI_FX / 2) + CAL_ONE).abs() < CAL_ONE / 200); // sin(3π/2) = −1
    }

    #[test]
    fn naive_pi_pulse_transfers_at_nominal() {
        let p = apply_composite(&[(0, PI_FX)], 0);
        assert!((p - CAL_ONE).abs() < CAL_ONE / 50, "naive P(0) = {p}");
    }

    #[test]
    fn robust_pulse_beats_naive() {
        let r = optimize_robust_pulse(fx(0.3), 13);
        // a naive π-pulse is fragile to 30% amplitude error
        assert!(r.worst_naive_fx > CAL_ONE / 12, "naive worst {} too low", r.worst_naive_fx);
        // the optimizer discovers a meaningfully more robust pulse
        assert!(
            r.worst_robust_fx * 2 < r.worst_naive_fx,
            "robust {} not better than naive {}",
            r.worst_robust_fx, r.worst_naive_fx
        );
        assert_eq!(r, optimize_robust_pulse(fx(0.3), 13)); // deterministic
    }

    #[cfg(feature = "quantum_ops")]
    #[test]
    fn robust_pulse_seals_a_receipt() {
        use crate::quantum_ops::GrantRef;
        use ed25519_dalek::SigningKey;
        let r = optimize_robust_pulse(fx(0.3), 13);
        let art = r.artifacts();
        let rec = seal_artifacts(
            &SigningKey::from_bytes(&[3u8; 32]), "did:key:lab", "sim:q0", &art,
            5_000_000, GrantRef::unbounded("quantum.calibrate"), None,
        );
        assert!(rec.verify());
        assert!(rec.evidence_matches(0, &art.evidence_bytes));
    }
}