prism-q 0.35.0

use super::*;
use crate::circuits;
use crate::sim::Engine;

type EngineFn = fn(&Circuit, &NoiseModel, usize, u64) -> Result<ShotsResult>;

fn assert_ghz_noise_spread(run: EngineFn, n: usize, min_each: usize) {
    let mut circuit = circuits::ghz_circuit(n);
    circuit.measure_all();

    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.01);
    let result = run(&circuit, &noise, 1000, 42).unwrap();

    assert_eq!(result.shots.len(), 1000);
    assert_eq!(result.shots[0].len(), n);

    let all_zero: Vec<bool> = vec![false; n];
    let all_one: Vec<bool> = vec![true; n];
    let num_00 = result.shots.iter().filter(|s| **s == all_zero).count();
    let num_11 = result.shots.iter().filter(|s| **s == all_one).count();

    assert!(
        1000 - num_00 - num_11 > 0,
        "noise should produce non-GHZ outcomes"
    );
    assert!(
        num_00 > min_each,
        "should still have many |00...0> outcomes"
    );
    assert!(
        num_11 > min_each,
        "should still have many |11...1> outcomes"
    );
}

fn assert_ghz_zero_noise_coherent(run: EngineFn, n: usize) {
    let mut circuit = circuits::ghz_circuit(n);
    circuit.measure_all();

    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.0);
    let result = run(&circuit, &noise, 100, 42).unwrap();
    assert_eq!(
        result.coherent_fraction(),
        1.0,
        "GHZ with zero noise must be all-0 or all-1"
    );
}

fn assert_clifford_noise_varies(run: EngineFn, n: usize) {
    let mut circuit = circuits::clifford_heavy_circuit(n, 10, 42);
    circuit.measure_all();

    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.01);
    let result = run(&circuit, &noise, 100, 42).unwrap();

    assert_eq!(result.shots.len(), 100);
    assert_eq!(result.shots[0].len(), n);

    let unique: std::collections::HashSet<Vec<bool>> = result.shots.iter().cloned().collect();
    assert!(unique.len() > 1, "noise should produce varied outcomes");
}

#[test]
fn noisy_ghz_produces_varied_outcomes() {
    assert_ghz_noise_spread(run_shots_noisy, 10, 100);
}

#[test]
fn zero_noise_matches_noiseless() {
    assert_ghz_zero_noise_coherent(run_shots_noisy, 5);
}

#[test]
fn noise_model_length_matches_circuit() {
    let n = 10;
    let mut circuit = circuits::ghz_circuit(n);
    circuit.measure_all();

    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.001);
    assert_eq!(noise.after_gate.len(), circuit.instructions.len());
}

#[test]
fn compiled_noisy_stats_match_brute_force() {
    let n = 10;
    let mut circuit = circuits::ghz_circuit(n);
    circuit.measure_all();

    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.01);
    let num_shots = 10000;

    let brute = run_shots_noisy_brute_with(
        |s| Box::new(StabilizerBackend::new(s)),
        &circuit,
        &noise,
        num_shots,
        42,
    )
    .unwrap();
    let compiled = run_shots_noisy_compiled(&circuit, &noise, num_shots, 42).unwrap();

    let brute_frac = brute.coherent_fraction();
    let compiled_frac = compiled.coherent_fraction();

    assert!(
        (brute_frac - compiled_frac).abs() < 0.05,
        "coherent fraction should be similar: brute={brute_frac:.3}, compiled={compiled_frac:.3}"
    );
    assert!(
        brute_frac < 1.0 && compiled_frac < 1.0,
        "both should produce errors"
    );
}

#[test]
fn pauli_engines_share_observable_statistics() {
    let n = 6;
    let mut circuit = Circuit::new(n, n);
    circuit.add_gate(Gate::X, &[0]);
    circuit.add_gate(Gate::H, &[2]);
    circuit.add_gate(Gate::Cx, &[2, 3]);
    circuit.add_gate(Gate::Cx, &[3, 4]);
    circuit.add_gate(Gate::H, &[5]);
    circuit.measure_all();

    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.02);
    let num_shots = 10_000;
    let seed = 42;

    let analytic =
        crate::sim::homological::noisy_marginals_analytical(&circuit, &noise, seed).unwrap();

    let sv_factory = |s: u64| -> Box<dyn Backend> {
        Box::new(crate::backend::statevector::StatevectorBackend::new(s))
    };
    let engines: Vec<(&str, ShotsResult)> = vec![
        (
            "compiled",
            run_shots_noisy_compiled(&circuit, &noise, num_shots, seed).unwrap(),
        ),
        (
            "frame",
            run_shots_noisy_frame(&circuit, &noise, num_shots, seed).unwrap(),
        ),
        (
            "brute",
            run_shots_noisy_brute_with(
                |s| Box::new(StabilizerBackend::new(s)),
                &circuit,
                &noise,
                num_shots,
                seed,
            )
            .unwrap(),
        ),
        (
            "trajectory",
            crate::sim::trajectory::run_trajectories(
                sv_factory,
                &circuit,
                &noise,
                num_shots,
                seed,
                false,
                crate::sim::ResolvedBackend::Statevector,
            )
            .unwrap(),
        ),
    ];

    // 5 sigma at 10k shots: sigma <= 0.5 / sqrt(N) = 0.005.
    let marginal_tol = 0.025;
    for (name, result) in &engines {
        assert_eq!(result.shots.len(), num_shots);
        for (bit, &expected) in analytic.iter().enumerate() {
            let p = result.marginal(bit);
            assert!(
                (p - expected).abs() < marginal_tol,
                "{name} bit {bit}: marginal {p:.4} vs analytic {expected:.4}"
            );
        }
    }

    // The density-matrix backend carries no sampling noise, so it sits far
    // tighter to the analytic marginals than the sampled engines. The
    // residual is a model gap, not statistics: the analytic reference is the
    // Pauli-frame independent-error model (X and Y flips are independent
    // Bernoulli events), while the density matrix evolves the exact
    // depolarizing channel (mutually exclusive Paulis). The two agree to
    // first order in p and differ at O(p^2); at p = 0.02 that gap is < 3e-4.
    let dm_marginals = dm_noisy_marginals(&circuit, &noise, seed).unwrap();
    let dm_tol = 3e-4;
    for (bit, &expected) in analytic.iter().enumerate() {
        assert!(
            (dm_marginals[bit] - expected).abs() < dm_tol,
            "density_matrix bit {bit}: marginal {:.12} vs analytic {expected:.12}",
            dm_marginals[bit]
        );
    }

    let ghz_coherent = |result: &ShotsResult| -> f64 {
        result
            .shots
            .iter()
            .filter(|s| s[2] == s[3] && s[3] == s[4])
            .count() as f64
            / num_shots as f64
    };
    let correlator_tol = 0.04;
    let reference = ghz_coherent(&engines[0].1);
    for (name, result) in &engines[1..] {
        let value = ghz_coherent(result);
        assert!(
            (value - reference).abs() < correlator_tol,
            "{name} GHZ coherence {value:.4} vs compiled {reference:.4}"
        );
    }
}

#[test]
fn dm_expectation_matches_run_expectation_values() {
    use crate::{PauliAxis, PauliTerm};
    let mut circuit = Circuit::new(3, 0);
    circuit.add_gate(Gate::H, &[0]);
    circuit.add_gate(Gate::Ry(0.7), &[1]);
    circuit.add_gate(Gate::Cx, &[0, 1]);
    circuit.add_gate(Gate::Rx(1.1), &[2]);
    circuit.add_gate(Gate::Cz, &[1, 2]);
    circuit.add_gate(Gate::T, &[0]);

    let observables = vec![
        vec![PauliTerm::new(0, PauliAxis::Z)],
        vec![PauliTerm::new(2, PauliAxis::X)],
        vec![
            PauliTerm::new(0, PauliAxis::Z),
            PauliTerm::new(1, PauliAxis::Z),
        ],
        vec![
            PauliTerm::new(0, PauliAxis::Y),
            PauliTerm::new(2, PauliAxis::X),
        ],
    ];

    let reference = crate::sim::run_expectation_values(&circuit, &observables, 42).unwrap();
    let dm = density_matrix_expectation_values(&circuit, &observables, None, 42).unwrap();
    for (i, (a, b)) in dm.iter().zip(&reference).enumerate() {
        assert!(
            (a - b).abs() < 1e-12,
            "observable {i}: dm {a:.15} vs statevector {b:.15}"
        );
    }
}

#[test]
fn trajectory_expectations_converge_to_density_matrix_within_5_sigma() {
    use crate::backend::statevector::StatevectorBackend;
    use crate::sim::trajectory::run_trajectory_shot;
    use crate::{PauliAxis, PauliTerm};

    // Six-qubit seeded layered circuit with two noise layers.
    let mut circuit = Circuit::new(6, 0);
    for q in 0..6 {
        circuit.add_gate(Gate::H, &[q]);
    }
    for q in 0..5 {
        circuit.add_gate(Gate::Cx, &[q, q + 1]);
    }
    let mut rng = ChaCha8Rng::seed_from_u64(0xDEAD_BEEF);
    for q in 0..6 {
        let theta: f64 = rand::RngExt::random::<f64>(&mut rng) * std::f64::consts::TAU;
        circuit.add_gate(
            if q % 2 == 0 {
                Gate::Rz(theta)
            } else {
                Gate::Rx(theta)
            },
            &[q],
        );
    }

    // Noise after layer 1 (the first H) and layer 2 (first CX).
    let mut noise = NoiseModel {
        after_gate: vec![Vec::new(); circuit.instructions.len()],
        readout: vec![None; circuit.num_classical_bits],
    };
    for q in 0..6 {
        noise.after_gate[q].push(NoiseEvent {
            channel: NoiseChannel::Depolarizing { p: 0.02 },
            qubits: smallvec![q],
        });
    }
    for q in 0..6 {
        noise.after_gate[6 + q.min(4)].push(NoiseEvent {
            channel: NoiseChannel::AmplitudeDamping { gamma: 0.03 },
            qubits: smallvec![q],
        });
    }

    let observables = vec![
        vec![
            PauliTerm::new(0, PauliAxis::Z),
            PauliTerm::new(1, PauliAxis::Z),
        ],
        vec![
            PauliTerm::new(2, PauliAxis::X),
            PauliTerm::new(3, PauliAxis::X),
        ],
        vec![
            PauliTerm::new(4, PauliAxis::Z),
            PauliTerm::new(5, PauliAxis::Z),
        ],
    ];
    let weights = [0.5f64, 0.3, 0.2];

    let mu_dm =
        density_matrix_expectation_values(&circuit, &observables, Some(&noise), 42).unwrap();

    // Trajectory estimates: each shot is a pure state; the observable value
    // on that state is a sample whose mean converges to Tr(rho P).
    let masks: Vec<_> = observables
        .iter()
        .map(|obs| crate::sim::pauli_masks(obs, circuit.num_qubits).unwrap())
        .collect();
    let shots = 20_000usize;
    let mut sums = vec![0.0f64; observables.len()];
    let mut sq_sums = vec![0.0f64; observables.len()];
    let mut backend = StatevectorBackend::new(42);
    let readout = crate::sim::trajectory::written_readout(&circuit, &noise.readout);
    let thermal_rates = crate::sim::trajectory::prepare_thermal_rates(&noise);
    for i in 0..shots {
        let mut shot_rng = ChaCha8Rng::seed_from_u64(42u64.wrapping_add(i as u64));
        run_trajectory_shot(
            &mut backend,
            &circuit,
            &noise,
            &readout,
            &thermal_rates,
            &mut shot_rng,
        )
        .unwrap();
        let state = backend.state_vector();
        let norm: f64 = state.iter().map(|a| a.norm_sqr()).sum();
        for (k, &(x, z, y)) in masks.iter().enumerate() {
            let v = crate::sim::pauli_expectation_from_masks(state, x, z, y, norm);
            sums[k] += v;
            sq_sums[k] += v * v;
        }
    }

    for (k, &w) in weights.iter().enumerate() {
        let mean = sums[k] / shots as f64;
        let var = (sq_sums[k] / shots as f64 - mean * mean).max(0.0);
        let sigma = (var / shots as f64).sqrt();
        let diff = (mean - mu_dm[k]).abs();
        if sigma < 1e-12 {
            assert!(
                diff < 1e-10,
                "term {k} (w={w}): degenerate trajectory, mean {mean} vs dm {}",
                mu_dm[k]
            );
        } else {
            let z_score = diff / sigma;
            assert!(
                z_score < 5.0,
                "term {k} (w={w}): z={z_score:.2}, traj {mean:.6} vs dm {:.6}",
                mu_dm[k]
            );
        }
    }
}

#[test]
fn dm_named_channel_lowering_preserves_trace_and_decay() {
    use crate::backend::density_matrix::DensityMatrixBackend;

    let trace = |dm: &DensityMatrixBackend| -> f64 { dm.probabilities().unwrap().iter().sum() };
    let identity = [
        [Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)],
        [Complex64::new(0.0, 0.0), Complex64::new(1.0, 0.0)],
    ];
    let channels = [
        NoiseChannel::Pauli {
            px: 0.1,
            py: 0.05,
            pz: 0.2,
        },
        NoiseChannel::Depolarizing { p: 0.3 },
        NoiseChannel::AmplitudeDamping { gamma: 0.4 },
        NoiseChannel::PhaseDamping { gamma: 0.35 },
        NoiseChannel::ThermalRelaxation {
            t1: 50.0,
            t2: 40.0,
            gate_time: 10.0,
            excited_population: 0.0,
        },
        NoiseChannel::Custom {
            kraus: vec![identity],
        },
    ];
    for ch in &channels {
        let mut dm = DensityMatrixBackend::new(42);
        dm.init(1, 0).unwrap();
        dm.apply(&Instruction::Gate {
            gate: Gate::H,
            targets: smallvec![0],
        })
        .unwrap();
        dm.apply_1q_kraus(0, &kraus_1q(ch));
        assert!(
            (trace(&dm) - 1.0).abs() < 1e-12,
            "trace after {ch:?}: {}",
            trace(&dm)
        );
    }

    let (t1, t2, gt) = (50.0, 40.0, 10.0);
    let thermal = NoiseChannel::ThermalRelaxation {
        t1,
        t2,
        gate_time: gt,
        excited_population: 0.0,
    };
    let mut dm = DensityMatrixBackend::new(42);
    dm.init(1, 0).unwrap();
    dm.apply(&Instruction::Gate {
        gate: Gate::X,
        targets: smallvec![0],
    })
    .unwrap();
    dm.apply_1q_kraus(0, &kraus_1q(&thermal));
    let rho = dm.reduced_density_matrix_1q(0).unwrap();
    assert!(
        (rho[1][1].re - (-gt / t1).exp()).abs() < 1e-12,
        "T1 population decay: {rho:?}"
    );

    let mut dm = DensityMatrixBackend::new(42);
    dm.init(1, 0).unwrap();
    dm.apply(&Instruction::Gate {
        gate: Gate::H,
        targets: smallvec![0],
    })
    .unwrap();
    dm.apply_1q_kraus(0, &kraus_1q(&thermal));
    let rho = dm.reduced_density_matrix_1q(0).unwrap();
    assert!(
        (rho[0][1].norm() - 0.5 * (-gt / t2).exp()).abs() < 1e-12,
        "T2 coherence decay: {rho:?}"
    );

    // Amplitude and phase damping lowerings drive the exact analytic channel
    // through kraus_1q, not just trace preservation, on |+>.
    let prep_plus = |ch: &NoiseChannel| -> [[Complex64; 2]; 2] {
        let mut dm = DensityMatrixBackend::new(42);
        dm.init(1, 0).unwrap();
        dm.apply(&Instruction::Gate {
            gate: Gate::H,
            targets: smallvec![0],
        })
        .unwrap();
        dm.apply_1q_kraus(0, &kraus_1q(ch));
        dm.reduced_density_matrix_1q(0).unwrap()
    };

    let gamma = 0.4;
    let rho = prep_plus(&NoiseChannel::AmplitudeDamping { gamma });
    assert!(
        (rho[0][0].re - (0.5 + 0.5 * gamma)).abs() < 1e-12,
        "AD pop0: {rho:?}"
    );
    assert!(
        (rho[1][1].re - 0.5 * (1.0 - gamma)).abs() < 1e-12,
        "AD pop1: {rho:?}"
    );
    assert!(
        (rho[0][1].norm() - 0.5 * (1.0 - gamma).sqrt()).abs() < 1e-12,
        "AD coherence: {rho:?}"
    );

    let gamma = 0.35;
    let rho = prep_plus(&NoiseChannel::PhaseDamping { gamma });
    assert!(
        (rho[0][0].re - 0.5).abs() < 1e-12,
        "PD pop0 preserved: {rho:?}"
    );
    assert!(
        (rho[1][1].re - 0.5).abs() < 1e-12,
        "PD pop1 preserved: {rho:?}"
    );
    assert!(
        (rho[0][1].norm() - 0.5 * (1.0 - gamma).sqrt()).abs() < 1e-12,
        "PD coherence: {rho:?}"
    );
}

#[test]
fn compiled_noisy_clifford_produces_noise() {
    assert_clifford_noise_varies(run_shots_noisy, 20);
}

#[test]
fn compile_noisy_rejects_reset_circuits() {
    let mut circuit = Circuit::new(1, 1);
    circuit.add_reset(0);
    circuit.add_measure(0, 0);
    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.01);
    assert!(compile_noisy(&circuit, &noise, 42).is_err());
}

#[test]
fn compile_noisy_rejects_conditionals() {
    let mut circuit = Circuit::new(2, 2);
    circuit.add_gate(Gate::H, &[0]);
    circuit.add_measure(0, 0);
    circuit.instructions.push(Instruction::Conditional {
        condition: crate::circuit::ClassicalCondition::BitIsOne(0),
        gate: Gate::X,
        targets: crate::circuit::smallvec![1],
    });
    circuit.add_measure(1, 1);
    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.01);
    assert!(compile_noisy(&circuit, &noise, 42).is_err());
}

#[test]
fn run_shots_noisy_handles_reset_circuits() {
    let mut circuit = Circuit::new(1, 1);
    circuit.add_gate(Gate::X, &[0]);
    circuit.add_reset(0);
    circuit.add_measure(0, 0);
    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.0);
    let result = run_shots_noisy(&circuit, &noise, 32, 42).unwrap();
    assert!(result.shots.iter().all(|shot| !shot[0]));
}

#[test]
fn run_shots_noisy_handles_conditionals() {
    let mut circuit = Circuit::new(2, 2);
    circuit.add_gate(Gate::H, &[0]);
    circuit.add_measure(0, 0);
    circuit.instructions.push(Instruction::Conditional {
        condition: crate::circuit::ClassicalCondition::BitIsOne(0),
        gate: Gate::X,
        targets: crate::circuit::smallvec![1],
    });
    circuit.add_measure(1, 1);
    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.0);
    let result = run_shots_noisy(&circuit, &noise, 256, 42).unwrap();
    assert!(result.shots.iter().all(|shot| shot[0] == shot[1]));
}

#[test]
fn frame_ghz_100q_produces_varied_outcomes() {
    assert_ghz_noise_spread(run_shots_noisy_frame, 100, 50);
}

#[test]
fn frame_zero_noise_matches_noiseless_100q() {
    assert_ghz_zero_noise_coherent(run_shots_noisy_frame, 100);
}

#[test]
fn frame_stats_match_compiled_ghz() {
    let n = 100;
    let mut circuit = circuits::ghz_circuit(n);
    circuit.measure_all();

    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.01);
    let num_shots = 5000;

    let frame = run_shots_noisy_frame(&circuit, &noise, num_shots, 42).unwrap();
    let compiled = run_shots_noisy_compiled(&circuit, &noise, num_shots, 42).unwrap();

    let frame_coh = frame.coherent_fraction();
    let compiled_coh = compiled.coherent_fraction();

    assert!(
        (frame_coh - compiled_coh).abs() < 0.05,
        "coherent fraction should be similar: frame={frame_coh:.3}, compiled={compiled_coh:.3}"
    );
}

// Every qubit is read twice in one terminal block: q0 is random, q1 reads one,
// and the second read of q2 overwrites bit 0.
#[test]
fn frame_rereads_a_qubit_measured_twice() {
    let mut circuit = Circuit::new(3, 5);
    circuit.add_gate(Gate::H, &[0]);
    circuit.add_gate(Gate::X, &[1]);
    circuit.add_measure(0, 0);
    circuit.add_measure(1, 1);
    circuit.add_measure(2, 2);
    circuit.add_measure(0, 3);
    circuit.add_measure(1, 4);
    circuit.add_measure(2, 0);

    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.0);
    let frame = run_shots_noisy_frame(&circuit, &noise, 512, 42).unwrap();
    let mut ones = 0;
    for shot in &frame.shots {
        assert_eq!(shot[..3], [false, true, false], "{shot:?}");
        assert!(shot[4], "{shot:?}");
        ones += usize::from(shot[3]);
    }
    assert!(
        (128..384).contains(&ones),
        "q0 read one {ones} of 512 times"
    );
}

#[test]
fn frame_clifford_100q_produces_noise() {
    assert_clifford_noise_varies(run_shots_noisy_frame, 100);
}

#[test]
fn filtered_noisy_bell_pairs_matches_monolithic() {
    let n_pairs = 50;
    let mut circuit = circuits::independent_bell_pairs(n_pairs);
    circuit.measure_all();

    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.01);
    let seed = 42u64;

    let filtered =
        compile_noisy_filtered(&circuit, &noise, &circuit.independent_subsystems(), seed).unwrap();
    let monolithic = compile_noisy_monolithic(&circuit, &noise, seed).unwrap();

    assert_eq!(filtered.num_measurements, monolithic.num_measurements);
    assert_eq!(filtered.events.len(), monolithic.events.len());

    let mut filtered = filtered;
    let mut monolithic = monolithic;
    let num_shots = 10_000;
    let shots_f = filtered.sample_bulk(num_shots);
    let shots_m = monolithic.sample_bulk(num_shots);

    assert_eq!(shots_f.len(), num_shots);
    assert_eq!(shots_m.len(), num_shots);

    let mut agree_f = 0usize;
    let mut agree_m = 0usize;
    for shot in &shots_f {
        for pair in shot.chunks(2) {
            if pair[0] == pair[1] {
                agree_f += 1;
            }
        }
    }
    for shot in &shots_m {
        for pair in shot.chunks(2) {
            if pair[0] == pair[1] {
                agree_m += 1;
            }
        }
    }

    let total_pairs = num_shots * n_pairs;
    let agree_rate_f = agree_f as f64 / total_pairs as f64;
    let agree_rate_m = agree_m as f64 / total_pairs as f64;
    assert!(
        agree_rate_f > 0.95,
        "filtered agreement rate {agree_rate_f:.4} should be >0.95 with low noise"
    );
    assert!(
        agree_rate_m > 0.95,
        "monolithic agreement rate {agree_rate_m:.4} should be >0.95 with low noise"
    );
    assert!(
        (agree_rate_f - agree_rate_m).abs() < 0.02,
        "filtered ({agree_rate_f:.4}) and monolithic ({agree_rate_m:.4}) should have similar agreement rates"
    );
}

// The frame and compiled routes answer from the same distribution, so each is
// pinned off the result in the same test that reads the rates.
fn readout_route_circuit(pad: usize) -> Circuit {
    let mut circuit = Circuit::new(3, 3);
    circuit.add_gate(Gate::X, &[0]);
    circuit.add_gate(Gate::X, &[1]);
    for _ in 0..pad {
        circuit.add_gate(Gate::Z, &[2]);
    }
    circuit.add_measure(0, 0);
    circuit.add_measure(1, 1);
    circuit.add_measure(2, 2);
    circuit
}

fn bit_rate(shots: &[Vec<bool>], bit: usize, want: bool) -> f64 {
    shots.iter().filter(|s| s[bit] == want).count() as f64 / shots.len() as f64
}

// Rates far apart in both directions: applied before the reference bits are
// folded in, bit 0 would read p01 = 0.40 where p10 = 0.05 belongs. The bound is
// five sigma on a rate estimated from 20k draws, rounded up.
fn assert_asymmetric_readout_rates(circuit: &Circuit, engine: Engine) {
    let mut noise = NoiseModel::uniform_depolarizing(circuit, 0.0);
    noise.set_bit_readout_error(0, 0.40, 0.05);
    noise.set_bit_readout_error(2, 0.20, 0.60);

    let result = crate::sim::simulate(circuit)
        .noise(&noise)
        .seed(42)
        .shots(20_000)
        .unwrap();
    assert_eq!(result.metadata.engine, Some(engine));

    let flipped_one = bit_rate(&result.shots, 0, false);
    assert!(
        (flipped_one - 0.05).abs() < 0.02,
        "bit 0 measures 1 and takes p10 = 0.05, got {flipped_one}"
    );
    let flipped_zero = bit_rate(&result.shots, 2, true);
    assert!(
        (flipped_zero - 0.20).abs() < 0.02,
        "bit 2 measures 0 and takes p01 = 0.20, got {flipped_zero}"
    );
    assert!(
        result.shots.iter().all(|s| s[1]),
        "bit 1 carries no readout entry and must keep its noiseless value"
    );
}

#[test]
fn readout_rates_hold_on_the_frame_route() {
    assert_asymmetric_readout_rates(&readout_route_circuit(0), Engine::FrameSampler);
}

#[test]
fn readout_rates_hold_on_the_compiled_route() {
    assert_asymmetric_readout_rates(&readout_route_circuit(10), Engine::NoisyCompiledSampler);
}

fn entangled_route_circuit(pad: usize) -> Circuit {
    let mut circuit = Circuit::new(3, 3);
    circuit.add_gate(Gate::H, &[0]);
    circuit.add_gate(Gate::Cx, &[0, 1]);
    circuit.add_gate(Gate::Cx, &[1, 2]);
    for _ in 0..pad {
        circuit.add_gate(Gate::Z, &[2]);
    }
    circuit.add_measure(0, 0);
    circuit.add_measure(1, 1);
    circuit.add_measure(2, 2);
    circuit
}

// The trajectory engine draws readout per shot against the unpacked record,
// which shares no code with the packed walk. Per-bit sigma is at most 0.0035 at
// 20k draws and the difference of two such estimates at most 0.005, so 0.03 is
// six sigma on the difference.
fn assert_matches_trajectory_engine(circuit: &Circuit, engine: Engine) {
    let mut noise = NoiseModel::uniform_depolarizing(circuit, 0.02);
    noise.set_bit_readout_error(0, 0.10, 0.30);
    noise.set_bit_readout_error(1, 0.25, 0.05);
    noise.set_bit_readout_error(2, 0.40, 0.10);

    let clifford = crate::sim::simulate(circuit)
        .noise(&noise)
        .seed(42)
        .shots(20_000)
        .unwrap();
    assert_eq!(clifford.metadata.engine, Some(engine));
    let trajectory = crate::sim::simulate(circuit)
        .backend(BackendKind::Statevector)
        .noise(&noise)
        .seed(42)
        .shots(20_000)
        .unwrap();
    assert_eq!(
        trajectory.metadata.backend,
        crate::sim::ResolvedBackend::Statevector
    );

    for bit in 0..3 {
        let c = bit_rate(&clifford.shots, bit, true);
        let t = bit_rate(&trajectory.shots, bit, true);
        assert!(
            (c - t).abs() < 0.03,
            "bit {bit}: clifford {c:.4} vs trajectory {t:.4}"
        );
    }
}

#[test]
fn readout_matches_the_trajectory_engine_on_the_frame_route() {
    assert_matches_trajectory_engine(&entangled_route_circuit(0), Engine::FrameSampler);
}

#[test]
fn readout_matches_the_trajectory_engine_on_the_compiled_route() {
    assert_matches_trajectory_engine(&entangled_route_circuit(7), Engine::NoisyCompiledSampler);
}

#[cfg(feature = "gpu")]
#[test]
fn noisy_gpu_test_circuit_routes_to_gpu_bts() {
    let mut circuit = Circuit::new(32, 32);
    for q in 0..16 {
        circuit.add_gate(Gate::H, &[q]);
    }
    for q in 16..32 {
        for k in 0..4 {
            circuit.add_gate(Gate::Cx, &[(q - 16 + k) % 16, q]);
        }
    }
    circuit.measure_all();
    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.02);

    let gpu = compile_noisy(&circuit, &noise, 42)
        .unwrap()
        .with_gpu(crate::gpu::GpuContext::stub_for_tests());
    assert!(
        gpu.noiseless
            .should_use_gpu_bts(crate::gpu::bts_min_shots()),
        "the golden_gpu noisy reduction test relies on this circuit shape \
         passing the GPU BTS routing gates"
    );
}

#[cfg(feature = "gpu")]
#[test]
fn compiled_noisy_with_stub_gpu_matches_cpu_below_threshold() {
    let n = 12;
    let mut circuit = circuits::ghz_circuit(n);
    circuit.measure_all();

    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.01);
    let shots = 20_000;

    let mut cpu = compile_noisy(&circuit, &noise, 42).unwrap();
    let cpu_marginals = cpu.sample_marginals(shots);

    let mut gpu = compile_noisy(&circuit, &noise, 42)
        .unwrap()
        .with_gpu(crate::gpu::GpuContext::stub_for_tests());
    let gpu_marginals = gpu.sample_marginals(shots);

    for (idx, (cpu_p1, gpu_p1)) in cpu_marginals.iter().zip(gpu_marginals.iter()).enumerate() {
        assert!(
            (cpu_p1 - gpu_p1).abs() < 0.03,
            "marginal[{idx}] diverged too much: cpu={cpu_p1}, gpu={gpu_p1}"
        );
    }
}

#[cfg(feature = "gpu")]
#[test]
fn compiled_noisy_with_stub_gpu_low_rank_above_threshold_uses_cpu_fallback() {
    let shots = crate::gpu::bts_min_shots().max(1);
    let n = 12;
    let mut circuit = circuits::ghz_circuit(n);
    circuit.measure_all();

    let noise = NoiseModel::uniform_depolarizing(&circuit, 0.01);
    let mut gpu = compile_noisy(&circuit, &noise, 42)
        .unwrap()
        .with_gpu(crate::gpu::GpuContext::stub_for_tests());

    assert!(!gpu.noiseless.should_use_gpu_bts(shots));
    let counts = gpu.sample_counts(shots);
    let total: u64 = counts.values().sum();
    assert_eq!(total, shots as u64);
}

// The identity product is not a branch: an index of 0 would leave 1 firing
// shot in 15 unchanged and cut the effective rate by exactly that much.
#[test]
fn pair_branch_covers_the_fifteen_non_identity_products() {
    let mut seen = [0usize; 16];
    for i in 0..15_000 {
        let (l0, l1) = pair_branch(i as f64 / 15_000.0);
        seen[l0 * 4 + l1] += 1;
    }
    assert_eq!(seen[0], 0, "identity is not a branch");
    assert!(seen[1..].iter().all(|&count| count == 1000), "{seen:?}");
}

fn pair_event(qubits: [usize; 2], p: f64) -> NoiseEvent {
    NoiseEvent {
        channel: NoiseChannel::TwoQubitDepolarizing { p },
        qubits: qubits.into_iter().collect(),
    }
}

fn bare_model(circuit: &Circuit) -> NoiseModel {
    NoiseModel {
        after_gate: vec![Vec::new(); circuit.instructions.len()],
        readout: vec![None; circuit.num_classical_bits],
    }
}

/// `X q0` then a CX chain, so the ideal record is all ones and every
/// single-qubit Z moves with the noise. `pad` is a run of Z on qubit 0, which
/// no Z-basis outcome sees and which moves the gate count per qubit across the
/// frame cutoff.
fn pair_chain_circuit(n: usize, pad: usize) -> Circuit {
    let mut circuit = Circuit::new(n, n);
    circuit.add_gate(Gate::X, &[0]);
    for q in 0..n - 1 {
        circuit.add_gate(Gate::Cx, &[q, q + 1]);
    }
    for _ in 0..pad {
        circuit.add_gate(Gate::Z, &[0]);
    }
    for q in 0..n {
        circuit.add_measure(q, q);
    }
    circuit
}

fn pair_after_every_cx(circuit: &Circuit, p: f64) -> NoiseModel {
    let mut model = bare_model(circuit);
    for (slot, inst) in model.after_gate.iter_mut().zip(&circuit.instructions) {
        if let Instruction::Gate {
            gate: Gate::Cx,
            targets,
        } = inst
        {
            slot.push(pair_event([targets[0], targets[1]], p));
        }
    }
    model
}

fn chain_observables(n: usize) -> Vec<Vec<crate::PauliTerm>> {
    let mut observables: Vec<Vec<crate::PauliTerm>> =
        (0..n).map(|q| vec![crate::PauliTerm::z(q)]).collect();
    observables
        .extend((0..n - 1).map(|q| vec![crate::PauliTerm::z(q), crate::PauliTerm::z(q + 1)]));
    observables
}

fn sampled_chain_expectations(shots: &[Vec<bool>], n: usize) -> Vec<f64> {
    let total = shots.len() as f64;
    let mut got: Vec<f64> = (0..n)
        .map(|q| 1.0 - 2.0 * shots.iter().filter(|s| s[q]).count() as f64 / total)
        .collect();
    got.extend(
        (0..n - 1)
            .map(|q| 1.0 - 2.0 * shots.iter().filter(|s| s[q] != s[q + 1]).count() as f64 / total),
    );
    got
}

/// Five sigma on a Z expectation estimated from 20k draws at the worst-case
/// variance, rounded up.
const CHAIN_BAND: f64 = 0.036;

/// `circuit` and `noise` without the terminal measurements and their event
/// slots, the unitary half the density-matrix oracle accepts.
fn unmeasured(circuit: &Circuit, noise: &NoiseModel) -> (Circuit, NoiseModel) {
    let after_gate = circuit
        .instructions
        .iter()
        .zip(&noise.after_gate)
        .filter(|(inst, _)| !matches!(inst, Instruction::Measure { .. }))
        .map(|(_, events)| events.clone())
        .collect();
    let noise = NoiseModel {
        after_gate,
        readout: noise.readout.clone(),
    };
    (circuit.without_measurements(), noise)
}

fn assert_chain_matches_density_matrix(
    circuit: &Circuit,
    noise: &NoiseModel,
    n: usize,
    engine: Engine,
) {
    let result = crate::sim::simulate(circuit)
        .noise(noise)
        .seed(42)
        .shots(20_000)
        .unwrap();
    assert_eq!(result.metadata.engine, Some(engine));
    let (unitary, unitary_noise) = unmeasured(circuit, noise);
    let want = density_matrix_expectation_values(
        &unitary,
        &chain_observables(n),
        Some(&unitary_noise),
        42,
    )
    .unwrap();

    for (i, (g, w)) in sampled_chain_expectations(&result.shots, n)
        .iter()
        .zip(&want)
        .enumerate()
    {
        assert!(
            (g - w).abs() < CHAIN_BAND,
            "observable {i}: sampled {g:.4} against density matrix {w:.4}"
        );
    }
}

// Every qubit accumulates one letter per event that reaches it, so a single
// qubit Z here separates the joint draw from an implementation that resolves
// only the first target. The nearest-neighbour parities cannot: each is fed by
// one letter of each of two events, and dropping the second leaves the same
// marginal.
#[test]
fn pair_channel_matches_the_density_matrix_on_the_frame_route() {
    let circuit = pair_chain_circuit(8, 0);
    assert_chain_matches_density_matrix(
        &circuit,
        &pair_after_every_cx(&circuit, 0.05),
        8,
        Engine::FrameSampler,
    );
}

#[test]
fn pair_channel_matches_the_density_matrix_on_the_compiled_route() {
    let circuit = pair_chain_circuit(8, 18);
    assert_chain_matches_density_matrix(
        &circuit,
        &pair_after_every_cx(&circuit, 0.05),
        8,
        Engine::NoisyCompiledSampler,
    );
}

// The grouped route XORs the single-qubit rows through a flip lookup table
// while the pair table runs as a second pass over the same buffer, so a model
// carrying both is the only thing that exercises the two together.
#[test]
fn a_mixed_model_runs_the_flip_table_and_the_pair_pass() {
    let n = 9;
    let circuit = pair_chain_circuit(n, 18);

    let mut noise = pair_after_every_cx(&circuit, 0.02);
    for (slot, inst) in noise.after_gate.iter_mut().zip(&circuit.instructions) {
        if let Instruction::Gate {
            gate: Gate::X | Gate::Cx,
            targets,
        } = inst
        {
            slot.extend(targets.iter().map(|&q| NoiseEvent {
                channel: NoiseChannel::Depolarizing { p: 0.055 },
                qubits: smallvec![q],
            }));
        }
    }

    let sampler = compile_noisy(&circuit, &noise, 42).unwrap();
    assert!(sampler.z_lut.is_some(), "the flip table must be built");
    assert!(sampler.events.has_pairs());

    assert_chain_matches_density_matrix(&circuit, &noise, n, Engine::NoisyCompiledSampler);
}

/// `|++>` under one two-qubit channel, read back in the X basis. `pad` is an
/// even run of H on qubit 0, which leaves the state alone and moves the gate
/// count per qubit across the frame cutoff.
fn pair_x_basis_circuit(pad: usize) -> Circuit {
    let mut circuit = Circuit::new(2, 2);
    for _ in 0..pad {
        circuit.add_gate(Gate::H, &[0]);
    }
    circuit.add_gate(Gate::H, &[0]);
    circuit.add_gate(Gate::H, &[1]);
    circuit.add_gate(Gate::H, &[0]);
    circuit.add_gate(Gate::H, &[1]);
    circuit.add_measure(0, 0);
    circuit.add_measure(1, 1);
    circuit
}

fn x_basis_pair_model(circuit: &Circuit, pad: usize, p: f64) -> NoiseModel {
    let mut model = bare_model(circuit);
    model.after_gate[pad + 1] = vec![pair_event([0, 1], p)];
    model
}

// In the X basis a bit reads 1 under Y or Z, so the outcome separates the 15
// branches into 3 that leave both bits at 0 and 4 that set both. Running the
// branch index 0 through 15 instead of 1 through 15 swaps identity in for ZZ
// and swaps those two rates. No Z-basis statistic can see that swap, because
// identity and ZZ flip the same measurements, namely none.
fn assert_identity_is_not_a_branch(circuit: &Circuit, pad: usize, engine: Engine) {
    let noise = x_basis_pair_model(circuit, pad, 1.0);
    let result = crate::sim::simulate(circuit)
        .noise(&noise)
        .seed(42)
        .shots(20_000)
        .unwrap();
    assert_eq!(result.metadata.engine, Some(engine));

    let total = result.shots.len() as f64;
    let both_zero = result.shots.iter().filter(|s| !s[0] && !s[1]).count() as f64 / total;
    let both_one = result.shots.iter().filter(|s| s[0] && s[1]).count() as f64 / total;
    assert!(
        (both_zero - 3.0 / 15.0).abs() < 0.02,
        "3 of the 15 branches leave both bits at 0, got {both_zero}"
    );
    assert!(
        (both_one - 4.0 / 15.0).abs() < 0.02,
        "4 of the 15 branches set both bits, got {both_one}"
    );
}

#[test]
fn identity_is_not_a_branch_on_the_frame_route() {
    assert_identity_is_not_a_branch(&pair_x_basis_circuit(0), 0, Engine::FrameSampler);
}

#[test]
fn identity_is_not_a_branch_on_the_compiled_route() {
    assert_identity_is_not_a_branch(&pair_x_basis_circuit(4), 4, Engine::NoisyCompiledSampler);
}

// One branch of 15 decides both letters, so each of the three outcomes that
// moves at least one bit carries `4p/15 = 0.133`. Two independent draws
// matching the per-qubit rate of `8p/15` would put `(8p/15)^2 = 0.071` on the
// corner where both bits move, which is the whole difference between a
// correlated channel and two uncorrelated ones.
fn assert_the_two_letters_are_drawn_jointly(circuit: &Circuit, engine: Engine) {
    let p = 0.5;
    let mut noise = bare_model(circuit);
    noise.after_gate[1] = vec![pair_event([0, 1], p)];

    let result = crate::sim::simulate(circuit)
        .noise(&noise)
        .seed(42)
        .shots(20_000)
        .unwrap();
    assert_eq!(result.metadata.engine, Some(engine));
    let total = result.shots.len() as f64;
    let rate = |want: [bool; 2]| {
        result
            .shots
            .iter()
            .filter(|s| s[0] == want[0] && s[1] == want[1])
            .count() as f64
            / total
    };

    let corner = 4.0 * p / 15.0;
    for (label, want, expected) in [
        ("both bits move", [false, false], corner),
        ("only qubit 0 moves", [false, true], corner),
        ("only qubit 1 moves", [true, false], corner),
        ("neither moves", [true, true], 1.0 - 3.0 * corner),
    ] {
        let got = rate(want);
        assert!(
            (got - expected).abs() < 0.02,
            "{label}: got {got:.4}, expected {expected:.4}"
        );
    }
}

#[test]
fn the_two_letters_are_drawn_jointly_on_the_frame_route() {
    assert_the_two_letters_are_drawn_jointly(&pair_chain_circuit(2, 0), Engine::FrameSampler);
}

#[test]
fn the_two_letters_are_drawn_jointly_on_the_compiled_route() {
    assert_the_two_letters_are_drawn_jointly(
        &pair_chain_circuit(2, 6),
        Engine::NoisyCompiledSampler,
    );
}

// A zero-rate pair flips nothing, so it must not cost the model the device
// path the way a live one does, and must not block the chain complex, which
// otherwise reads the channel variant alone. This is the reading the readout
// table already uses: inert, not absent.
#[test]
fn a_zero_rate_pair_channel_is_inert() {
    let circuit = pair_x_basis_circuit(4);

    let inert = x_basis_pair_model(&circuit, 4, 0.0);
    let sampler = compile_noisy(&circuit, &inert, 42).unwrap();
    assert!(!sampler.events.has_pairs());
    assert!(inert.is_pauli_only());
    assert!(crate::sim::homological::noisy_marginals_analytical(&circuit, &inert, 42).is_ok());

    let live = x_basis_pair_model(&circuit, 4, 0.2);
    let sampler = compile_noisy(&circuit, &live, 42).unwrap();
    assert!(sampler.events.has_pairs());
    assert!(!live.is_pauli_only());
}

// The inert entry must also draw nothing. A pass that draws and never fires
// leaves the record alone, so the draw shows only in what the stream feeds
// next: the readout entry, drawn from the same generator after the pair pass.
// The two runs share one seed and diverge on the first extra draw.
#[test]
fn a_zero_rate_pair_channel_consumes_no_randomness() {
    let mut circuit = Circuit::new(2, 2);
    circuit.add_gate(Gate::H, &[0]);
    circuit.add_gate(Gate::Cx, &[0, 1]);
    for _ in 0..6 {
        circuit.add_gate(Gate::Z, &[0]);
    }
    circuit.add_measure(0, 0);
    circuit.add_measure(1, 1);

    let mut singles = NoiseModel::uniform_depolarizing(&circuit, 0.08);
    singles.set_bit_readout_error(0, 0.3, 0.3);
    let mut with_inert_pair = NoiseModel::uniform_depolarizing(&circuit, 0.08);
    with_inert_pair.set_bit_readout_error(0, 0.3, 0.3);
    with_inert_pair.after_gate[1].push(pair_event([0, 1], 0.0));

    let baseline = run_shots_noisy(&circuit, &singles, 999, 42).unwrap();
    let inert = run_shots_noisy(&circuit, &with_inert_pair, 999, 42).unwrap();
    assert_eq!(baseline.metadata.engine, Some(Engine::NoisyCompiledSampler));
    assert_eq!(inert.metadata.engine, Some(Engine::NoisyCompiledSampler));
    assert!(
        baseline.shots.iter().any(|s| s[0]) && baseline.shots.iter().any(|s| !s[0]),
        "the fixture must carry a random outcome for the comparison to mean anything"
    );
    assert_eq!(baseline.shots, inert.shots);
}

/// Two subsystems that no gate couples, over a deterministic all-ones record.
fn two_block_circuit() -> Circuit {
    let mut circuit = Circuit::new(4, 4);
    circuit.add_gate(Gate::X, &[0]);
    circuit.add_gate(Gate::Cx, &[0, 1]);
    circuit.add_gate(Gate::X, &[2]);
    circuit.add_gate(Gate::Cx, &[2, 3]);
    for q in 0..4 {
        circuit.add_measure(q, q);
    }
    circuit
}

// A pair whose two qubits land in different subsystem blocks has no single
// block to be propagated in, so the model takes the monolithic compile. A pair
// inside one block keeps the filtered one.
#[test]
fn a_block_straddling_pair_leaves_the_filtered_compile() {
    let circuit = two_block_circuit();
    let blocks = circuit.independent_subsystems();
    assert_eq!(blocks.len(), 2);

    let model = |qubits: [usize; 2], p: f64| {
        let mut noise = bare_model(&circuit);
        noise.after_gate[1] = vec![pair_event(qubits, p)];
        noise
    };

    assert!(pairs_stay_in_one_block(
        &circuit,
        &model([0, 1], 0.2),
        &blocks
    ));
    assert!(!pairs_stay_in_one_block(
        &circuit,
        &model([1, 2], 0.2),
        &blocks
    ));
    assert!(
        pairs_stay_in_one_block(&circuit, &model([1, 2], 0.0), &blocks),
        "a zero-rate pair stores no row and constrains no route"
    );
}

// The filtered compile scatters each block's local masks back into the global
// measurement order, four rows for a pair where a single-qubit event has two.
// Nothing else reaches that scatter carrying a live pair.
#[test]
fn the_filtered_compile_carries_an_in_block_pair() {
    let circuit = two_block_circuit();
    let blocks = circuit.independent_subsystems();
    let mut noise = bare_model(&circuit);
    noise.after_gate[1] = vec![pair_event([0, 1], 0.5)];
    assert!(pairs_stay_in_one_block(&circuit, &noise, &blocks));

    let mut sampler = compile_noisy_filtered(&circuit, &noise, &blocks, 42).unwrap();
    assert!(sampler.events.has_pairs());
    let shots = sampler.sample_bulk_packed(20_000).to_shots();

    let observables: Vec<Vec<crate::PauliTerm>> = (0..4)
        .map(|q| vec![crate::PauliTerm::z(q)])
        .chain(std::iter::once(vec![
            crate::PauliTerm::z(0),
            crate::PauliTerm::z(1),
        ]))
        .collect();
    let (unitary, unitary_noise) = unmeasured(&circuit, &noise);
    let want = density_matrix_expectation_values(&unitary, &observables, Some(&unitary_noise), 42)
        .unwrap();

    let total = shots.len() as f64;
    let mut got: Vec<f64> = (0..4)
        .map(|q| 1.0 - 2.0 * shots.iter().filter(|s| s[q]).count() as f64 / total)
        .collect();
    got.push(1.0 - 2.0 * shots.iter().filter(|s| s[0] != s[1]).count() as f64 / total);

    for (i, (g, w)) in got.iter().zip(&want).enumerate() {
        assert!(
            (g - w).abs() < CHAIN_BAND,
            "observable {i}: sampled {g:.4} against density matrix {w:.4}"
        );
    }
}

// The exact walk defers one-qubit maps and folds them, a gate, and its
// channels into one 16x16 sweep. Applying every instruction and channel on its
// own is the reference: orientation (a CX with its control above its target),
// non-unital damping, noise on one target only, a gate with nothing to fold, a
// channel with no gate, and a three-qubit gate that flushes all run through it.
#[test]
fn deferred_density_matrix_walk_matches_one_sweep_per_map() {
    use crate::backend::density_matrix::DensityMatrixBackend;

    // Four qubits take the native walk, seven fold selectively, nine fold any
    // pending map.
    for n in [4, 7, 9] {
        let mut circuit = Circuit::new(n, 0);
        let mut after_gate: Vec<Vec<NoiseEvent>> = Vec::new();
        let damp = |q: usize| NoiseEvent {
            channel: NoiseChannel::AmplitudeDamping { gamma: 0.07 },
            qubits: smallvec::smallvec![q],
        };
        let mut push = |gate: Gate, targets: &[usize], events: Vec<NoiseEvent>| {
            circuit.add_gate(gate, targets);
            after_gate.push(events);
        };
        push(Gate::H, &[0], vec![NoiseEvent::pauli(0, 0.01, 0.02, 0.03)]);
        push(Gate::Ry(0.4), &[1], vec![damp(1)]);
        push(Gate::Rx(1.1), &[1], vec![]);
        push(
            Gate::Cx,
            &[1, 0],
            vec![NoiseEvent::pauli(1, 0.02, 0.0, 0.01), damp(0)],
        );
        push(Gate::T, &[2], vec![]);
        push(Gate::Rzz(0.8), &[2, 3], vec![damp(3)]);
        push(Gate::Swap, &[0, 2], vec![]);
        push(Gate::Cz, &[3, 1], vec![]);
        push(
            Gate::cu(Gate::Ry(0.9).matrix_2x2()),
            &[0, 3],
            vec![damp(0), damp(3)],
        );
        push(Gate::S, &[1], vec![NoiseEvent::pauli(2, 0.05, 0.0, 0.0)]);
        push(
            Gate::mcu(Gate::X.matrix_2x2(), 2),
            &[0, 1, 2],
            vec![damp(2)],
        );
        push(
            Gate::Rz(0.3),
            &[3],
            vec![NoiseEvent::pauli(3, 0.0, 0.04, 0.0)],
        );
        let noise = NoiseModel {
            after_gate,
            readout: Vec::new(),
        };

        let folded = evolve_density_matrix(
            &BackendKind::DensityMatrix,
            &circuit,
            Some(&noise),
            None,
            42,
        )
        .unwrap()
        .density_matrix()
        .unwrap();

        let mut reference = DensityMatrixBackend::new(42);
        reference.init(n, 0).unwrap();
        for (inst, events) in circuit.instructions.iter().zip(&noise.after_gate) {
            reference.apply(inst).unwrap();
            for event in events {
                apply_noise_event_dm(&mut reference, event);
            }
        }
        let reference = reference.density_matrix().unwrap();

        for (i, (a, b)) in folded.iter().zip(&reference).enumerate() {
            assert!(
                (a - b).norm() < 1e-12,
                "{n} qubits, entry {i}: folded {a} vs reference {b}"
            );
        }
    }
}

/// Apply every instruction and channel on its own, measurements skipped as the
/// exact walk skips them.
fn native_density_matrix_walk(
    circuit: &Circuit,
    noise: &NoiseModel,
) -> crate::backend::density_matrix::DensityMatrixBackend {
    let mut reference = crate::backend::density_matrix::DensityMatrixBackend::new(42);
    reference
        .init(circuit.num_qubits, circuit.num_classical_bits)
        .unwrap();
    for (inst, events) in circuit.instructions.iter().zip(&noise.after_gate) {
        if !matches!(inst, Instruction::Measure { .. }) {
            reference.apply(inst).unwrap();
        }
        for event in events {
            apply_noise_event_dm(&mut reference, event);
        }
    }
    reference
}

/// The folded walk against [`native_density_matrix_walk`]: every matrix entry,
/// then one-qubit X, Y, and Z expectations and an X(x)Y pair on each neighbour.
fn assert_walk_matches_native(circuit: &Circuit, noise: &NoiseModel, label: &str) {
    use crate::PauliTerm;

    let n = circuit.num_qubits;
    let folded =
        evolve_density_matrix(&BackendKind::DensityMatrix, circuit, Some(noise), None, 42).unwrap();
    let reference = native_density_matrix_walk(circuit, noise);
    let (a, b) = (
        folded.density_matrix().unwrap(),
        reference.density_matrix().unwrap(),
    );
    for (i, (x, y)) in a.iter().zip(&b).enumerate() {
        assert!(
            (x - y).norm() < 1e-12,
            "{label}, entry {i}: folded {x} vs reference {y}"
        );
    }

    let mut observables: Vec<Vec<PauliTerm>> = (0..n)
        .flat_map(|q| {
            [
                vec![PauliTerm::x(q)],
                vec![PauliTerm::y(q)],
                vec![PauliTerm::z(q)],
            ]
        })
        .collect();
    observables.extend((0..n - 1).map(|q| vec![PauliTerm::x(q), PauliTerm::y(q + 1)]));
    let masks: Vec<_> = observables
        .iter()
        .map(|obs| crate::sim::pauli_masks(obs, n).unwrap())
        .collect();
    let (x, y) = (
        folded.expectations_pauli(&masks),
        reference.expectations_pauli(&masks),
    );
    for (i, (x, y)) in x.iter().zip(&y).enumerate() {
        assert!(
            (x - y).abs() < 1e-12,
            "{label}, observable {i}: folded {x} vs reference {y}"
        );
    }
}

// Two-qubit channels on a gate's own pair fold into its sweep in either
// orientation, one-qubit channels on other qubits wait on those qubits, and a
// channel straddling the pair ends the fold. Barriers carry idle channels and
// keep pending maps alive.
#[test]
fn exact_walk_folds_pair_and_idle_channels() {
    let damp = |q: usize, gamma: f64| NoiseEvent {
        channel: NoiseChannel::AmplitudeDamping { gamma },
        qubits: smallvec::smallvec![q],
    };
    let relax = |q: usize| NoiseEvent {
        channel: NoiseChannel::ThermalRelaxation {
            t1: 50.0,
            t2: 40.0,
            gate_time: 1.5,
            excited_population: 0.0,
        },
        qubits: smallvec::smallvec![q],
    };
    let depol2 = |a: usize, b: usize, p: f64| NoiseEvent {
        channel: NoiseChannel::TwoQubitDepolarizing { p },
        qubits: smallvec::smallvec![a, b],
    };
    // {sqrt(1-p) I, sqrt(p) (X(x)Z) with a phase on one entry}: not symmetric
    // under swapping the pair, so a wrong orientation shows.
    let kraus2 = |a: usize, b: usize| {
        let z = Complex64::new(0.0, 0.0);
        let p: f64 = 0.06;
        let mut k0 = [[z; 4]; 4];
        let mut k1 = [[z; 4]; 4];
        for t in 0..4 {
            k0[t][t] = Complex64::new((1.0 - p).sqrt(), 0.0);
            let sign = if t & 1 == 0 { 1.0 } else { -1.0 };
            k1[t ^ 2][t] = Complex64::new(p.sqrt() * sign, 0.0);
        }
        k1[3][1] *= Complex64::new(0.0, 1.0);
        NoiseEvent {
            channel: NoiseChannel::Kraus2q {
                kraus: vec![k0, k1],
            },
            qubits: smallvec::smallvec![a, b],
        }
    };

    for n in [4, 7, 9] {
        let mut circuit = Circuit::new(n, 0);
        let mut after_gate: Vec<Vec<NoiseEvent>> = Vec::new();
        let mut push = |gate: Gate, targets: &[usize], events: Vec<NoiseEvent>| {
            circuit.add_gate(gate, targets);
            after_gate.push(events);
        };
        for q in 0..4 {
            push(Gate::Ry(0.3 + 0.4 * q as f64), &[q], vec![]);
        }
        push(Gate::H, &[0], vec![damp(2, 0.05), relax(3)]);
        push(
            Gate::Cx,
            &[0, 1],
            vec![relax(0), relax(1), depol2(0, 1, 0.04), damp(3, 0.02)],
        );
        push(Gate::Cz, &[2, 1], vec![kraus2(1, 2), damp(2, 0.1)]);
        push(Gate::Cx, &[3, 2], vec![kraus2(3, 2), relax(0)]);
        push(
            Gate::Fused2q(Box::new(Gate::Rzz(0.7).matrix_4x4())),
            &[1, 3],
            vec![
                damp(1, 0.03),
                kraus2(3, 1),
                relax(3),
                depol2(1, 3, 0.02),
                damp(0, 0.04),
            ],
        );
        push(
            Gate::cu(Gate::Ry(0.9).matrix_2x2()),
            &[0, 2],
            vec![
                damp(0, 0.05),
                kraus2(2, 3),
                damp(0, 0.07),
                depol2(0, 2, 0.03),
            ],
        );
        push(Gate::Rzz(0.5), &[1, 2], vec![depol2(2, 1, 0.05)]);
        push(Gate::Swap, &[0, 3], vec![]);
        push(Gate::Sdg, &[1], vec![depol2(0, 2, 0.01), damp(1, 0.02)]);
        circuit.add_barrier(&(0..n).collect::<Vec<_>>());
        after_gate.push(vec![relax(0), relax(2), depol2(1, 3, 0.02)]);
        let mut push = |gate: Gate, targets: &[usize], events: Vec<NoiseEvent>| {
            circuit.add_gate(gate, targets);
            after_gate.push(events);
        };
        push(Gate::Rx(0.8), &[2], vec![relax(1)]);
        push(
            Gate::Cx,
            &[2, 0],
            vec![relax(2), relax(0), depol2(2, 0, 0.03)],
        );
        push(
            Gate::mcu(Gate::X.matrix_2x2(), 2),
            &[0, 1, 3],
            vec![damp(2, 0.04)],
        );
        push(Gate::Cz, &[1, 2], vec![kraus2(1, 2), relax(3)]);
        let noise = NoiseModel {
            after_gate,
            readout: Vec::new(),
        };
        assert_walk_matches_native(&circuit, &noise, &format!("{n} qubits"));
    }
}

/// `n` qubits of `H`, `T`, `Rz` per layer, then a CX chain, measured.
fn non_clifford_layers(n: usize, depth: usize) -> Circuit {
    let mut circuit = Circuit::new(n, n);
    for layer in 0..depth {
        for q in 0..n {
            circuit.add_gate(Gate::H, &[q]);
            circuit.add_gate(Gate::T, &[q]);
            circuit.add_gate(Gate::Rz(0.03 * (layer + q + 1) as f64), &[q]);
        }
        for q in 0..n - 1 {
            circuit.add_gate(Gate::Cx, &[q, q + 1]);
        }
    }
    for q in 0..n {
        circuit.add_measure(q, q);
    }
    circuit
}

#[test]
fn exact_walk_matches_native_on_device_presets() {
    use crate::sim::calibration::presets;

    for n in [4, 7, 8] {
        let circuit = non_clifford_layers(n, 2);
        for (name, calibration) in [
            ("transmon", presets::superconducting_transmon(n)),
            ("neutral atom", presets::neutral_atom(n)),
        ] {
            let plain = calibration.to_noise_model(&circuit).unwrap();
            assert!(plain.readout.iter().any(Option::is_some));
            assert_walk_matches_native(&circuit, &plain, &format!("{name}, {n} qubits"));
            let scheduled = calibration.to_scheduled_noise_model(&circuit).unwrap();
            assert_walk_matches_native(
                &circuit,
                &scheduled,
                &format!("{name} scheduled, {n} qubits"),
            );
        }
    }
}

#[test]
fn noiseless_exact_walk_runs_the_fused_stream() {
    use crate::PauliTerm;

    for n in [5, 8] {
        let circuit = circuits::random_circuit(n, 6, 42);
        let observables: Vec<Vec<PauliTerm>> = (0..n)
            .flat_map(|q| [vec![PauliTerm::x(q)], vec![PauliTerm::y(q)]])
            .collect();
        let fused = density_matrix_expectation_values(&circuit, &observables, None, 42).unwrap();
        let mut reference = crate::backend::density_matrix::DensityMatrixBackend::new(42);
        reference.init(n, 0).unwrap();
        reference.apply_instructions(&circuit.instructions).unwrap();
        let masks: Vec<_> = observables
            .iter()
            .map(|obs| crate::sim::pauli_masks(obs, n).unwrap())
            .collect();
        for (i, (x, y)) in fused
            .iter()
            .zip(reference.expectations_pauli(&masks))
            .enumerate()
        {
            assert!(
                (x - y).abs() < 1e-12,
                "{n} qubits, observable {i}: fused {x} vs reference {y}"
            );
        }
    }
}

/// Each qubit's detuning over one unit-length layer, per shot, as the
/// trajectory engine draws it.
fn drawn_detunings(drift: DriftDistribution, num_qubits: usize, shots: usize) -> Vec<Vec<f64>> {
    let mut circuit = Circuit::new(num_qubits, 0);
    for q in 0..num_qubits {
        circuit.add_gate(Gate::Id, &[q]);
    }
    let noise = NoiseBuilder::new()
        .schedule(crate::sim::calibration::GateTimes::new(1.0, 1.0))
        .quasi_static_detuning(drift)
        .build(&circuit)
        .unwrap();
    let mut weights = vec![Vec::new(); num_qubits];
    for event in noise.after_gate.iter().flatten() {
        if let NoiseChannel::QuasiStatic { axis, weights: w } = &event.channel {
            assert_eq!(*axis, PauliAxis::Z);
            weights[event.qubits[0]] = w.clone();
        }
    }
    let mut z = vec![0.0; noise.quasi_static_sources()];
    (0..shots)
        .map(|shot| {
            let mut rng = crate::sim::trajectory::noise_rng(crate::sim::mix_seed(42, shot));
            crate::sim::trajectory::standard_normals(&mut rng, &mut z);
            weights
                .iter()
                .map(|w| w.iter().map(|&(source, weight)| weight * z[source]).sum())
                .collect()
        })
        .collect()
}

fn assert_covariance(samples: &[Vec<f64>], expected: &[f64]) {
    let n = samples[0].len();
    let count = samples.len() as f64;
    for a in 0..n {
        for b in 0..n {
            let observed = samples.iter().map(|s| s[a] * s[b]).sum::<f64>() / count;
            let want = expected[a * n + b];
            let spread = ((expected[a * n + a] * expected[b * n + b] + want * want) / count).sqrt();
            assert!(
                (observed - want).abs() <= 5.0 * spread + 1e-12,
                "covariance ({a}, {b}): drew {observed}, asked {want} +- {spread}"
            );
        }
    }
}

#[test]
fn correlated_detuning_draws_reproduce_the_requested_covariance() {
    let neighbours = DriftDistribution::independent([1.0, 2.0, 0.5, 1.5])
        .with_neighbour_correlation([(0, 1), (1, 2), (2, 3)], 0.45);
    let expected = neighbours.covariance().to_vec();
    assert_covariance(&drawn_detunings(neighbours, 4, 40_000), &expected);

    let full = DriftDistribution::from_covariance(vec![
        vec![2.0, 0.6, -0.4],
        vec![0.6, 1.0, 0.3],
        vec![-0.4, 0.3, 0.5],
    ]);
    let expected = full.covariance().to_vec();
    assert_covariance(&drawn_detunings(full, 3, 40_000), &expected);

    let locked =
        DriftDistribution::independent([1.0, 1.0]).with_neighbour_correlation([(0, 1)], 1.0);
    for shot in drawn_detunings(locked, 2, 64) {
        assert!((shot[0] - shot[1]).abs() < 1e-12, "{shot:?}");
    }
}