wai-quantum 0.3.26

//! Exact minimum-weight perfect-matching decoding — `wai.quantum.match`.
//!
//! A graphlike detector error model is a graph. Each detector is a node; one
//! more node is the boundary. Each error mechanism is an edge between the
//! detectors it flips, or between one detector and the boundary. An edge
//! weighs `ln((1 − p)/p)`: the less likely the error, the heavier the edge.
//!
//! The fired detectors of a shot must be explained by a set of errors. The
//! most likely set is a minimum-weight perfect matching:
//! - each fired detector pairs with another, or with the boundary;
//! - a pair costs the shortest path between its two ends;
//! - the predicted observable flip is the parity the matched paths carry.
//!
//! This module finds that matching **exactly**.
//! - Shortest paths come from Dijkstra on integer edge weights, so path costs
//!   are exact sums.
//! - The matching itself is Edmonds' weighted blossom algorithm, maximising
//!   `C − cost` over maximum-cardinality matchings. Fired detectors get
//!   boundary twins joined at zero cost, so a perfect matching always exists
//!   when the boundary is reachable.
//!
//! Two cuts keep the searches local without giving up exactness. Every
//! detector's path to the boundary is computed once per graph. A pair costing
//! at least both ends' boundary paths is never needed, because the zero-cost
//! twins absorb the swap; so each search stops at that cap.
//!
//! Everything is integer and ordered, so a decode is the same on every machine.
//! [`MatchingGraph::from_model`] refuses a mechanism with a part touching more
//! than two detectors, rather than dropping it: such a model is not graphlike,
//! and a matching decoder would silently mis-score it.
//!
//! # What it measures
//!
//! [`memory_experiment`] samples and decodes the rotated surface-code memory
//! experiment (`d` rounds, every noise knob at `p`). Logical error per round:
//!
//! | d | p = 0.002 | 0.003 | 0.004 | 0.005 |
//! |---|-----------|-------|-------|-------|
//! | 3 | 1.0e-3 | 2.2e-3 | 3.7e-3 | 5.6e-3 |
//! | 5 | 2.1e-4 | 6.6e-4 | 1.5e-3 | 2.8e-3 |
//! | 7 | 4.0e-5 | 1.9e-4 | 6.6e-4 | 1.4e-3 |
//!
//! Distance pays at every rate here, and [`fit_code_model`] turns the grid into
//! the model resource estimates assume, `P(d) = a·(p/p*)^((d+1)/2)` per round:
//! `a ≈ 0.023`, `p* ≈ 0.0100`. Every point was checked against an independent
//! sampler and decoder on the same circuits, and none differs by more than 1.4
//! standard errors.
//!
//! # Cost
//!
//! The decoder is exact first and fast second. At small distances it decodes
//! in microseconds. At `d = 9` a fired detector's boundary cap still spans
//! much of the lattice, so a shot costs about half a millisecond. A
//! region-growing matcher is the step after this one. Sampling is not the
//! bottleneck: the frame simulator draws over a million `d = 9` shots per second.

use crate::quantum_frame::{Detections, ErrorModel};
use crate::repro::ln;
use std::cmp::Reverse;
use std::collections::BinaryHeap;

/// Integer units per unit of log-likelihood weight.
const SCALE: f64 = 16_384.0;

#[derive(Clone, Debug, PartialEq, Eq)]
pub enum MatchError {
    /// A mechanism part touches this many detectors (more than two).
    Hyperedge(usize),
    /// No perfect matching: a fired detector cannot reach a partner or the
    /// boundary.
    Unmatchable,
}

/// The matching graph of a detector error model.
#[derive(Clone, Debug, PartialEq)]
pub struct MatchingGraph {
    detectors: usize,
    /// Neighbours of each node (detectors, then the boundary at index
    /// `detectors`): `(node, integer weight, observable mask)`.
    adj: Vec<Vec<(u32, i64, u64)>>,
    /// Probability mass of parts that flip observables but no detector:
    /// logical errors no decoder can see.
    pub undetectable: f64,
    /// Each detector's shortest path to the boundary: (cost, observable
    /// parity), or `None` if it cannot reach it.
    to_boundary: Vec<Option<(i64, u64)>>,
}

impl MatchingGraph {
    /// Build the graph from the model's per-component parts. Parallel edges
    /// with the same observable effect are XOR-combined; otherwise the more
    /// likely one is kept.
    pub fn from_model(m: &ErrorModel) -> Result<MatchingGraph, MatchError> {
        let n = m.detectors as usize;
        let boundary = n as u32;
        // (u, v) -> (p, obs), u < v.
        let mut edges: std::collections::BTreeMap<(u32, u32), (f64, u64)> = std::collections::BTreeMap::new();
        let mut undetectable = 0.0;
        for mech in &m.mechanisms {
            let p = mech.probability;
            for (dets, obs) in &mech.parts {
                let key = match dets.as_slice() {
                    [] => {
                        if *obs != 0 {
                            undetectable = undetectable + p - 2.0 * undetectable * p;
                        }
                        continue;
                    }
                    [a] => (*a, boundary),
                    [a, b] => (*a.min(b), *a.max(b)),
                    more => return Err(MatchError::Hyperedge(more.len())),
                };
                edges
                    .entry(key)
                    .and_modify(|e| {
                        if e.1 == *obs {
                            e.0 = e.0 + p - 2.0 * e.0 * p;
                        } else if p > e.0 {
                            *e = (p, *obs);
                        }
                    })
                    .or_insert((p, *obs));
            }
        }
        let mut adj = vec![Vec::new(); n + 1];
        for ((u, v), (p, obs)) in edges {
            let w = if p >= 0.5 { 0 } else { (ln((1.0 - p) / p) * SCALE).round() as i64 };
            adj[u as usize].push((v, w, obs));
            adj[v as usize].push((u, w, obs));
        }
        let mut g = MatchingGraph { detectors: n, adj, undetectable, to_boundary: Vec::new() };
        g.to_boundary = g.paths_from(boundary, None)[..n].to_vec();
        Ok(g)
    }

    /// Shortest paths from `src`, as (cost, observable parity), never passing
    /// through the boundary. With `targets` (each with a cost cap), the search
    /// stops once every target is settled or no unsettled target can still be
    /// reached within its cap.
    fn paths_from(&self, src: u32, targets: Option<&[(u32, i64)]>) -> Vec<Option<(i64, u64)>> {
        let n = self.detectors + 1;
        let boundary = self.detectors as u32;
        let mut best: Vec<Option<(i64, u64)>> = vec![None; n];
        let mut done = vec![false; n];
        let mut open: Vec<(u32, i64)> = targets.map_or(Vec::new(), |t| t.to_vec());
        let mut limit = open.iter().map(|t| t.1).max().unwrap_or(i64::MAX);
        let mut heap = BinaryHeap::new();
        best[src as usize] = Some((0, 0));
        heap.push(Reverse((0i64, src, 0u64)));
        while let Some(Reverse((c, u, obs))) = heap.pop() {
            if done[u as usize] {
                continue;
            }
            if targets.is_some() && c >= limit {
                break;
            }
            done[u as usize] = true;
            if let Some(i) = open.iter().position(|t| t.0 == u) {
                open.swap_remove(i);
                if open.is_empty() {
                    break;
                }
                limit = open.iter().map(|t| t.1).max().unwrap();
            }
            if u == boundary && u != src {
                continue;
            }
            for &(v, w, o) in &self.adj[u as usize] {
                // Per-shot searches leave the boundary to the precomputed table.
                if targets.is_some() && v == boundary {
                    continue;
                }
                let nc = c + w;
                if !done[v as usize] && best[v as usize].is_none_or(|(bc, _)| nc < bc) {
                    best[v as usize] = Some((nc, obs ^ o));
                    heap.push(Reverse((nc, v, obs ^ o)));
                }
            }
        }
        best
    }

    /// Predict which observables flipped, from the detectors that fired.
    pub fn decode(&self, fired: &[u32]) -> Result<u64, MatchError> {
        let k = fired.len();
        if k == 0 {
            return Ok(0);
        }
        let reach = |i: usize| self.to_boundary[fired[i] as usize];
        let mut edges: Vec<(usize, usize, i64)> = Vec::new();
        let mut parity: Vec<u64> = Vec::new();
        let mut max_cost = 0;
        for i in 0..k {
            // A pair costing at least both ends' boundary paths never helps:
            // the zero-cost boundary twins absorb the swap. So the search from
            // `i` looks for each later `j` only within that cap.
            let cap = |j: usize| match (reach(i), reach(j)) {
                (Some((a, _)), Some((b, _))) => a + b,
                _ => i64::MAX,
            };
            let targets: Vec<(u32, i64)> = (i + 1..k).map(|j| (fired[j], cap(j))).collect();
            let paths = if targets.is_empty() { Vec::new() } else { self.paths_from(fired[i], Some(&targets)) };
            for j in i + 1..k {
                if let Some((c, o)) = paths[fired[j] as usize] {
                    let dominated = match (reach(i), reach(j)) {
                        (Some((a, _)), Some((b, _))) => c >= a + b,
                        _ => false,
                    };
                    if !dominated {
                        edges.push((i, j, c));
                        parity.push(o);
                        max_cost = max_cost.max(c);
                    }
                }
            }
            if let Some((c, o)) = reach(i) {
                edges.push((i, k + i, c));
                parity.push(o);
                max_cost = max_cost.max(c);
            }
            for j in i + 1..k {
                edges.push((k + i, k + j, 0));
                parity.push(0);
            }
        }
        // Maximise C − cost over maximum-cardinality matchings.
        let big = max_cost + 1;
        let weighted: Vec<(usize, usize, i64)> = edges.iter().map(|&(a, b, c)| (a, b, big - c)).collect();
        let mate = max_weight_matching(2 * k, &weighted, true);
        if mate.iter().any(|m| m.is_none()) {
            return Err(MatchError::Unmatchable);
        }
        let mut out = 0;
        for (e, &(a, b, _)) in edges.iter().enumerate() {
            if mate[a] == Some(b) {
                out ^= parity[e];
            }
        }
        Ok(out)
    }

    /// Decode every shot; returns how many predictions missed the actual
    /// observable flips.
    pub fn failures(&self, shots: &Detections) -> Result<u64, MatchError> {
        let mut n = 0;
        for s in 0..shots.shots {
            if self.decode(&shots.fired(s))? != shots.flips(s) {
                n += 1;
            }
        }
        Ok(n)
    }
}

// ---------------------------------------------------------------------------
// Memory experiments and the code model they measure
// ---------------------------------------------------------------------------

/// The outcome of one surface-code memory experiment.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct MemoryPoint {
    pub distance: u32,
    pub rounds: u32,
    /// Physical error rate (every noise knob).
    pub p: f64,
    pub shots: u64,
    /// Shots whose decoded observable was wrong.
    pub failures: u64,
}

impl MemoryPoint {
    /// Logical error per round, `(1 − (1 − 2P)^{1/r}) / 2` for a per-shot
    /// failure rate `P` over `r` rounds.
    pub fn per_round(&self) -> f64 {
        let shot = (self.failures as f64 / self.shots as f64).min(0.499_999);
        (1.0 - crate::repro::exp(ln(1.0 - 2.0 * shot) / self.rounds as f64)) / 2.0
    }
}

/// Sample and decode `shots` runs of the rotated surface-code memory
/// experiment (distance `d`, `rounds` rounds, every noise knob at `p`).
pub fn memory_experiment(d: u32, rounds: u32, p: f64, shots: usize, seed: u64) -> Result<MemoryPoint, MatchError> {
    use crate::quantum_frame::{error_model, parse, sample, surface_code_memory};
    let c = parse(&surface_code_memory(d, rounds, p)).expect("the generated circuit parses");
    let g = MatchingGraph::from_model(&error_model(&c).expect("the generated circuit is deterministic"))?;
    let failures = g.failures(&sample(&c, shots, seed))?;
    Ok(MemoryPoint { distance: d, rounds, p, shots: shots as u64, failures })
}

/// A code model `P(d) = a·(p/p*)^((d+1)/2)` per round, as measured.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct CodeFit {
    /// The prefactor `a`.
    pub prefactor: f64,
    /// The threshold `p*`.
    pub threshold: f64,
}

#[cfg(feature = "quantum_resource")]
impl CodeFit {
    /// The gate-based surface-code model of a resource estimate, with this
    /// fit's constants in place of assumed ones.
    pub fn surface_code(&self) -> crate::quantum_resource::QecModel {
        crate::quantum_resource::QecModel {
            name: "surface code (measured)".into(),
            prefactor: self.prefactor,
            threshold: self.threshold,
            ..crate::quantum_resource::QecModel::surface_gate()
        }
    }
}

/// Fit the code model to measured points by least squares on
/// `ln P − x·ln p = ln a − x·ln p*`, with `x = (d+1)/2`. Each point is
/// weighted by its failure count, the inverse variance of `ln P` to first
/// order. Points without failures carry no information and are skipped.
/// `None` without two distinct distances.
pub fn fit_code_model(points: &[MemoryPoint]) -> Option<CodeFit> {
    let pts: Vec<(f64, f64, f64)> = points
        .iter()
        .filter(|q| q.failures > 0)
        .map(|q| {
            let x = f64::from(q.distance + 1) / 2.0;
            (x, ln(q.per_round()) - x * ln(q.p), q.failures as f64)
        })
        .collect();
    let wsum: f64 = pts.iter().map(|t| t.2).sum();
    if pts.is_empty() || wsum == 0.0 {
        return None;
    }
    let xm = pts.iter().map(|t| t.2 * t.0).sum::<f64>() / wsum;
    let ym = pts.iter().map(|t| t.2 * t.1).sum::<f64>() / wsum;
    let sxx: f64 = pts.iter().map(|t| t.2 * (t.0 - xm) * (t.0 - xm)).sum();
    if sxx == 0.0 {
        return None;
    }
    let sxy: f64 = pts.iter().map(|t| t.2 * (t.0 - xm) * (t.1 - ym)).sum();
    let slope = sxy / sxx;
    let intercept = ym - slope * xm;
    Some(CodeFit { prefactor: crate::repro::exp(intercept), threshold: crate::repro::exp(-slope) })
}

// ---------------------------------------------------------------------------
// Edmonds' weighted blossom algorithm, O(n³)
// ---------------------------------------------------------------------------

const NONE: usize = usize::MAX;

/// A maximum-weight matching of an undirected graph on `n` vertices with
/// integer edge weights. With `max_cardinality`, the maximum weight among
/// maximum-cardinality matchings. Returns each vertex's mate.
pub fn max_weight_matching(n: usize, edges: &[(usize, usize, i64)], max_cardinality: bool) -> Vec<Option<usize>> {
    if edges.is_empty() || n == 0 {
        return vec![None; n];
    }
    Blossom::new(n, edges).solve(max_cardinality)
}

struct Blossom<'a> {
    n: usize,
    edges: &'a [(usize, usize, i64)],
    endpoint: Vec<usize>,
    neighbend: Vec<Vec<usize>>,
    mate: Vec<usize>,
    label: Vec<u8>,
    labelend: Vec<usize>,
    inblossom: Vec<usize>,
    blossomparent: Vec<usize>,
    blossomchilds: Vec<Vec<usize>>,
    blossombase: Vec<usize>,
    blossomendps: Vec<Vec<usize>>,
    bestedge: Vec<usize>,
    blossombestedges: Vec<Option<Vec<usize>>>,
    unusedblossoms: Vec<usize>,
    dualvar: Vec<i64>,
    allowedge: Vec<bool>,
    queue: Vec<usize>,
}

impl<'a> Blossom<'a> {
    fn new(n: usize, edges: &'a [(usize, usize, i64)]) -> Self {
        let ne = edges.len();
        let maxweight = edges.iter().map(|e| e.2).max().unwrap_or(0).max(0);
        let mut endpoint = Vec::with_capacity(2 * ne);
        let mut neighbend = vec![Vec::new(); n];
        for (k, &(i, j, _)) in edges.iter().enumerate() {
            endpoint.push(i);
            endpoint.push(j);
            neighbend[i].push(2 * k + 1);
            neighbend[j].push(2 * k);
        }
        let mut dualvar = vec![maxweight; n];
        dualvar.extend(std::iter::repeat_n(0, n));
        Blossom {
            n,
            edges,
            endpoint,
            neighbend,
            mate: vec![NONE; n],
            label: vec![0; 2 * n],
            labelend: vec![NONE; 2 * n],
            inblossom: (0..n).collect(),
            blossomparent: vec![NONE; 2 * n],
            blossomchilds: vec![Vec::new(); 2 * n],
            blossombase: (0..n).chain(std::iter::repeat_n(NONE, n)).collect(),
            blossomendps: vec![Vec::new(); 2 * n],
            bestedge: vec![NONE; 2 * n],
            blossombestedges: vec![None; 2 * n],
            unusedblossoms: (n..2 * n).collect(),
            dualvar,
            allowedge: vec![false; ne],
            queue: Vec::new(),
        }
    }

    fn slack(&self, k: usize) -> i64 {
        let (i, j, w) = self.edges[k];
        self.dualvar[i] + self.dualvar[j] - 2 * w
    }

    fn leaves(&self, b: usize, out: &mut Vec<usize>) {
        if b < self.n {
            out.push(b);
        } else {
            for &t in &self.blossomchilds[b] {
                self.leaves(t, out);
            }
        }
    }

    fn leaves_of(&self, b: usize) -> Vec<usize> {
        let mut v = Vec::new();
        self.leaves(b, &mut v);
        v
    }

    fn assign_label(&mut self, w: usize, t: u8, p: usize) {
        let b = self.inblossom[w];
        self.label[w] = t;
        self.label[b] = t;
        self.labelend[w] = p;
        self.labelend[b] = p;
        self.bestedge[w] = NONE;
        self.bestedge[b] = NONE;
        if t == 1 {
            let l = self.leaves_of(b);
            self.queue.extend(l);
        } else if t == 2 {
            let base = self.blossombase[b];
            let mb = self.mate[base];
            self.assign_label(self.endpoint[mb], 1, mb ^ 1);
        }
    }

    fn scan_blossom(&mut self, mut v: usize, mut w: usize) -> usize {
        let mut path = Vec::new();
        let mut base = NONE;
        while v != NONE || w != NONE {
            let mut b = self.inblossom[v];
            if self.label[b] & 4 != 0 {
                base = self.blossombase[b];
                break;
            }
            path.push(b);
            self.label[b] = 5;
            if self.labelend[b] == NONE {
                v = NONE;
            } else {
                v = self.endpoint[self.labelend[b]];
                b = self.inblossom[v];
                v = self.endpoint[self.labelend[b]];
            }
            if w != NONE {
                core::mem::swap(&mut v, &mut w);
            }
        }
        for b in path {
            self.label[b] = 1;
        }
        base
    }

    fn add_blossom(&mut self, base: usize, k: usize) {
        let (mut v, mut w, _) = self.edges[k];
        let bb = self.inblossom[base];
        let mut bv = self.inblossom[v];
        let mut bw = self.inblossom[w];
        let b = self.unusedblossoms.pop().expect("a free blossom slot");
        self.blossombase[b] = base;
        self.blossomparent[b] = NONE;
        self.blossomparent[bb] = b;
        let mut path = Vec::new();
        let mut endps = Vec::new();
        while bv != bb {
            self.blossomparent[bv] = b;
            path.push(bv);
            endps.push(self.labelend[bv]);
            v = self.endpoint[self.labelend[bv]];
            bv = self.inblossom[v];
        }
        path.push(bb);
        path.reverse();
        endps.reverse();
        endps.push(2 * k);
        while bw != bb {
            self.blossomparent[bw] = b;
            path.push(bw);
            endps.push(self.labelend[bw] ^ 1);
            w = self.endpoint[self.labelend[bw]];
            bw = self.inblossom[w];
        }
        self.label[b] = 1;
        self.labelend[b] = self.labelend[bb];
        self.dualvar[b] = 0;
        self.blossomchilds[b] = path.clone();
        self.blossomendps[b] = endps;
        for v in self.leaves_of(b) {
            if self.label[self.inblossom[v]] == 2 {
                self.queue.push(v);
            }
            self.inblossom[v] = b;
        }
        let mut bestedgeto = vec![NONE; 2 * self.n];
        for &bv in &path {
            let lists: Vec<Vec<usize>> = match self.blossombestedges[bv].take() {
                None => self.leaves_of(bv).iter().map(|&v| self.neighbend[v].iter().map(|p| p / 2).collect()).collect(),
                Some(l) => vec![l],
            };
            for list in lists {
                for k in list {
                    let (mut i, mut j, _) = self.edges[k];
                    if self.inblossom[j] == b {
                        core::mem::swap(&mut i, &mut j);
                    }
                    let _ = i;
                    let bj = self.inblossom[j];
                    if bj != b
                        && self.label[bj] == 1
                        && (bestedgeto[bj] == NONE || self.slack(k) < self.slack(bestedgeto[bj]))
                    {
                        bestedgeto[bj] = k;
                    }
                }
            }
            self.bestedge[bv] = NONE;
        }
        let best: Vec<usize> = bestedgeto.into_iter().filter(|&k| k != NONE).collect();
        self.bestedge[b] = NONE;
        for &k in &best {
            if self.bestedge[b] == NONE || self.slack(k) < self.slack(self.bestedge[b]) {
                self.bestedge[b] = k;
            }
        }
        self.blossombestedges[b] = Some(best);
    }

    fn expand_blossom(&mut self, b: usize, endstage: bool) {
        let childs = self.blossomchilds[b].clone();
        for &s in &childs {
            self.blossomparent[s] = NONE;
            if s < self.n {
                self.inblossom[s] = s;
            } else if endstage && self.dualvar[s] == 0 {
                self.expand_blossom(s, endstage);
            } else {
                for v in self.leaves_of(s) {
                    self.inblossom[v] = s;
                }
            }
        }
        if !endstage && self.label[b] == 2 {
            let entrychild = self.inblossom[self.endpoint[self.labelend[b] ^ 1]];
            let len = childs.len() as isize;
            let mut j = childs.iter().position(|&c| c == entrychild).unwrap() as isize;
            let (jstep, endptrick): (isize, usize) = if j & 1 == 1 {
                j -= len;
                (1, 0)
            } else {
                (-1, 1)
            };
            let at = |j: isize| -> usize { (j.rem_euclid(len)) as usize };
            let endps = self.blossomendps[b].clone();
            let mut p = self.labelend[b];
            while j != 0 {
                let ep = self.endpoint[p ^ 1];
                self.label[ep] = 0;
                let q = endps[at(j - endptrick as isize)] ^ endptrick ^ 1;
                self.label[self.endpoint[q]] = 0;
                self.assign_label(ep, 2, p);
                self.allowedge[endps[at(j - endptrick as isize)] / 2] = true;
                j += jstep;
                p = endps[at(j - endptrick as isize)] ^ endptrick;
                self.allowedge[p / 2] = true;
                j += jstep;
            }
            let bv = childs[at(j)];
            let ep = self.endpoint[p ^ 1];
            self.label[ep] = 2;
            self.label[bv] = 2;
            self.labelend[ep] = p;
            self.labelend[bv] = p;
            self.bestedge[bv] = NONE;
            j += jstep;
            while childs[at(j)] != entrychild {
                let bv = childs[at(j)];
                if self.label[bv] == 1 {
                    j += jstep;
                    continue;
                }
                let mut found = NONE;
                for v in self.leaves_of(bv) {
                    if self.label[v] != 0 {
                        found = v;
                        break;
                    }
                }
                if found != NONE {
                    let v = found;
                    self.label[v] = 0;
                    let mb = self.mate[self.blossombase[bv]];
                    self.label[self.endpoint[mb]] = 0;
                    let le = self.labelend[v];
                    self.assign_label(v, 2, le);
                }
                j += jstep;
            }
        }
        self.label[b] = 0;
        self.labelend[b] = NONE;
        self.blossomchilds[b].clear();
        self.blossomendps[b].clear();
        self.blossombase[b] = NONE;
        self.blossombestedges[b] = None;
        self.bestedge[b] = NONE;
        self.unusedblossoms.push(b);
    }

    fn augment_blossom(&mut self, b: usize, v: usize) {
        let mut t = v;
        while self.blossomparent[t] != b {
            t = self.blossomparent[t];
        }
        if t >= self.n {
            self.augment_blossom(t, v);
        }
        let childs = self.blossomchilds[b].clone();
        let endps = self.blossomendps[b].clone();
        let len = childs.len() as isize;
        let i = childs.iter().position(|&c| c == t).unwrap();
        let mut j = i as isize;
        let (jstep, endptrick): (isize, usize) = if i & 1 == 1 {
            j -= len;
            (1, 0)
        } else {
            (-1, 1)
        };
        let at = |j: isize| -> usize { (j.rem_euclid(len)) as usize };
        while j != 0 {
            j += jstep;
            let t = childs[at(j)];
            let p = endps[at(j - endptrick as isize)] ^ endptrick;
            if t >= self.n {
                self.augment_blossom(t, self.endpoint[p]);
            }
            j += jstep;
            let t = childs[at(j)];
            if t >= self.n {
                self.augment_blossom(t, self.endpoint[p ^ 1]);
            }
            self.mate[self.endpoint[p]] = p ^ 1;
            self.mate[self.endpoint[p ^ 1]] = p;
        }
        let mut c = childs[i..].to_vec();
        c.extend_from_slice(&childs[..i]);
        let mut e = endps[i..].to_vec();
        e.extend_from_slice(&endps[..i]);
        self.blossombase[b] = self.blossombase[c[0]];
        self.blossomchilds[b] = c;
        self.blossomendps[b] = e;
    }

    fn augment_matching(&mut self, k: usize) {
        let (v, w, _) = self.edges[k];
        for (mut s, mut p) in [(v, 2 * k + 1), (w, 2 * k)] {
            loop {
                let bs = self.inblossom[s];
                if bs >= self.n {
                    self.augment_blossom(bs, s);
                }
                self.mate[s] = p;
                if self.labelend[bs] == NONE {
                    break;
                }
                let t = self.endpoint[self.labelend[bs]];
                let bt = self.inblossom[t];
                s = self.endpoint[self.labelend[bt]];
                let j = self.endpoint[self.labelend[bt] ^ 1];
                if bt >= self.n {
                    self.augment_blossom(bt, j);
                }
                self.mate[j] = self.labelend[bt];
                p = self.labelend[bt] ^ 1;
            }
        }
    }

    fn solve(mut self, max_cardinality: bool) -> Vec<Option<usize>> {
        let n = self.n;
        for _ in 0..n {
            self.label.iter_mut().for_each(|l| *l = 0);
            self.bestedge.iter_mut().for_each(|e| *e = NONE);
            for b in n..2 * n {
                self.blossombestedges[b] = None;
            }
            self.allowedge.iter_mut().for_each(|a| *a = false);
            self.queue.clear();
            for v in 0..n {
                if self.mate[v] == NONE && self.label[self.inblossom[v]] == 0 {
                    self.assign_label(v, 1, NONE);
                }
            }
            let mut augmented = false;
            loop {
                while let Some(v) = self.queue.pop() {
                    if augmented {
                        break;
                    }
                    for idx in 0..self.neighbend[v].len() {
                        let p = self.neighbend[v][idx];
                        let k = p / 2;
                        let w = self.endpoint[p];
                        if self.inblossom[v] == self.inblossom[w] {
                            continue;
                        }
                        let mut kslack = 0;
                        if !self.allowedge[k] {
                            kslack = self.slack(k);
                            if kslack <= 0 {
                                self.allowedge[k] = true;
                            }
                        }
                        if self.allowedge[k] {
                            if self.label[self.inblossom[w]] == 0 {
                                self.assign_label(w, 2, p ^ 1);
                            } else if self.label[self.inblossom[w]] == 1 {
                                let base = self.scan_blossom(v, w);
                                if base != NONE {
                                    self.add_blossom(base, k);
                                } else {
                                    self.augment_matching(k);
                                    augmented = true;
                                    break;
                                }
                            } else if self.label[w] == 0 {
                                self.label[w] = 2;
                                self.labelend[w] = p ^ 1;
                            }
                        } else if self.label[self.inblossom[w]] == 1 {
                            let b = self.inblossom[v];
                            if self.bestedge[b] == NONE || kslack < self.slack(self.bestedge[b]) {
                                self.bestedge[b] = k;
                            }
                        } else if self.label[w] == 0
                            && (self.bestedge[w] == NONE || kslack < self.slack(self.bestedge[w]))
                        {
                            self.bestedge[w] = k;
                        }
                    }
                }
                if augmented {
                    break;
                }
                // Choose the dual adjustment.
                let mut deltatype = 0u8;
                let mut delta = 0i64;
                let mut deltaedge = NONE;
                let mut deltablossom = NONE;
                if !max_cardinality {
                    deltatype = 1;
                    delta = self.dualvar[..n].iter().copied().min().unwrap();
                }
                for v in 0..n {
                    if self.label[self.inblossom[v]] == 0 && self.bestedge[v] != NONE {
                        let d = self.slack(self.bestedge[v]);
                        if deltatype == 0 || d < delta {
                            delta = d;
                            deltatype = 2;
                            deltaedge = self.bestedge[v];
                        }
                    }
                }
                for b in 0..2 * n {
                    if self.blossomparent[b] == NONE && self.label[b] == 1 && self.bestedge[b] != NONE {
                        let kslack = self.slack(self.bestedge[b]);
                        debug_assert_eq!(kslack % 2, 0);
                        let d = kslack / 2;
                        if deltatype == 0 || d < delta {
                            delta = d;
                            deltatype = 3;
                            deltaedge = self.bestedge[b];
                        }
                    }
                }
                for b in n..2 * n {
                    if self.blossombase[b] != NONE
                        && self.blossomparent[b] == NONE
                        && self.label[b] == 2
                        && (deltatype == 0 || self.dualvar[b] < delta)
                    {
                        delta = self.dualvar[b];
                        deltatype = 4;
                        deltablossom = b;
                    }
                }
                if deltatype == 0 {
                    deltatype = 1;
                    delta = self.dualvar[..n].iter().copied().min().unwrap().max(0);
                }
                for v in 0..n {
                    match self.label[self.inblossom[v]] {
                        1 => self.dualvar[v] -= delta,
                        2 => self.dualvar[v] += delta,
                        _ => {}
                    }
                }
                for b in n..2 * n {
                    if self.blossombase[b] != NONE && self.blossomparent[b] == NONE {
                        match self.label[b] {
                            1 => self.dualvar[b] += delta,
                            2 => self.dualvar[b] -= delta,
                            _ => {}
                        }
                    }
                }
                match deltatype {
                    1 => break,
                    2 => {
                        self.allowedge[deltaedge] = true;
                        let (mut i, mut j, _) = self.edges[deltaedge];
                        if self.label[self.inblossom[i]] == 0 {
                            core::mem::swap(&mut i, &mut j);
                        }
                        let _ = j;
                        self.queue.push(i);
                    }
                    3 => {
                        self.allowedge[deltaedge] = true;
                        let (i, _, _) = self.edges[deltaedge];
                        self.queue.push(i);
                    }
                    _ => self.expand_blossom(deltablossom, false),
                }
            }
            if !augmented {
                break;
            }
            for b in n..2 * n {
                if self.blossomparent[b] == NONE
                    && self.blossombase[b] != NONE
                    && self.label[b] == 1
                    && self.dualvar[b] == 0
                {
                    self.expand_blossom(b, true);
                }
            }
        }
        self.mate.iter().map(|&p| if p == NONE { None } else { Some(self.endpoint[p]) }).collect()
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::quantum_frame::{error_model, parse};

    /// Minimum-cost perfect matching by exhaustion (None if there is none).
    fn brute(n: usize, cost: &[Vec<Option<i64>>]) -> Option<i64> {
        fn go(used: &mut Vec<bool>, cost: &[Vec<Option<i64>>]) -> Option<i64> {
            let Some(i) = used.iter().position(|u| !u) else { return Some(0) };
            used[i] = true;
            let mut best = None;
            for j in i + 1..used.len() {
                if used[j] {
                    continue;
                }
                if let Some(c) = cost[i][j] {
                    used[j] = true;
                    if let Some(rest) = go(used, cost) {
                        best = Some(best.map_or(c + rest, |b: i64| b.min(c + rest)));
                    }
                    used[j] = false;
                }
            }
            used[i] = false;
            best
        }
        go(&mut vec![false; n], cost)
    }

    #[test]
    #[allow(clippy::needless_range_loop)]
    fn blossom_finds_the_minimum_perfect_matching() {
        let mut state = 12345u64;
        let mut next = |m: u64| {
            state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
            (state >> 33) % m
        };
        for trial in 0..3000 {
            let n = 2 * (1 + next(5) as usize); // 2..10 vertices
            let density = 40 + next(61); // percent of pairs present
            let mut cost = vec![vec![None; n]; n];
            let mut edges = Vec::new();
            for i in 0..n {
                for j in i + 1..n {
                    if next(100) < density {
                        let c = next(30) as i64;
                        cost[i][j] = Some(c);
                        cost[j][i] = Some(c);
                        edges.push((i, j, c));
                    }
                }
            }
            let want = brute(n, &cost);
            let big = 1000;
            let w: Vec<(usize, usize, i64)> = edges.iter().map(|&(a, b, c)| (a, b, big - c)).collect();
            let mate = max_weight_matching(n, &w, true);
            match want {
                None => assert!(mate.iter().any(|m| m.is_none()), "trial {trial}: matched with no perfect matching"),
                Some(c) => {
                    assert!(mate.iter().all(|m| m.is_some()), "trial {trial}: not perfect");
                    let mut total = 0;
                    for (i, m) in mate.iter().enumerate() {
                        let j = m.unwrap();
                        assert_eq!(mate[j], Some(i));
                        if i < j {
                            total += cost[i][j].expect("a real edge");
                        }
                    }
                    assert_eq!(total, c, "trial {trial}");
                }
            }
        }
    }

    #[test]
    fn the_repetition_code_decodes_every_correctable_error() {
        // Distance 5, one round of perfect checks: every pattern of at most 2
        // data flips is corrected.
        let d = 5;
        let mut text = String::from("R 0 1 2 3 4 5 6 7 8\n");
        text += "X_ERROR(0.01) 0 1 2 3 4\n";
        for i in 0..d - 1 {
            text += &format!("CX {} {} {} {}\n", i, d + i, i + 1, d + i);
        }
        text += "MR 5 6 7 8\nM 0 1 2 3 4\n";
        for i in 0..d - 1 {
            text += &format!("DETECTOR rec[-{}]\n", 9 - i);
        }
        text += "OBSERVABLE_INCLUDE(0) rec[-1]\n";
        let c = parse(&text).unwrap();
        let g = MatchingGraph::from_model(&error_model(&c).unwrap()).unwrap();
        for mask in 0u32..32 {
            if mask.count_ones() > 2 {
                continue;
            }
            let flipped = |q: u32| mask >> q & 1;
            let fired: Vec<u32> = (0..d - 1).filter(|&i| flipped(i) ^ flipped(i + 1) == 1).collect();
            assert_eq!(g.decode(&fired).unwrap(), flipped(d - 1) as u64, "mask {mask:05b}");
        }
    }

    #[test]
    fn the_fit_recovers_a_known_model() {
        // Exact per-round rates from a = 0.07, p* = 0.009, turned into failure
        // counts over a huge number of shots (so rounding is negligible).
        let (a, ps): (f64, f64) = (0.07, 0.009);
        let mut pts = Vec::new();
        for d in [3u32, 5, 7] {
            for p in [0.002, 0.004] {
                let per_round: f64 = a * (p / ps).powi(d.div_ceil(2) as i32);
                let rounds = d;
                let shot = (1.0 - (1.0 - 2.0 * per_round).powi(rounds as i32)) / 2.0;
                let shots = 1u64 << 50;
                pts.push(MemoryPoint { distance: d, rounds, p, shots, failures: (shot * shots as f64).round() as u64 });
            }
        }
        let f = fit_code_model(&pts).unwrap();
        assert!((f.prefactor / a - 1.0).abs() < 1e-6, "{f:?}");
        assert!((f.threshold / ps - 1.0).abs() < 1e-6, "{f:?}");
        // One distance cannot separate a from p*.
        assert!(fit_code_model(&pts[..2]).is_none());
        #[cfg(feature = "quantum_resource")]
        {
            let m = f.surface_code();
            assert_eq!((m.prefactor, m.threshold), (f.prefactor, f.threshold));
            assert_eq!(m.tile_qubits(13), crate::quantum_resource::QecModel::surface_gate().tile_qubits(13));
        }
    }

    #[test]
    fn distance_pays_on_the_surface_code() {
        // Below threshold, d = 5 beats d = 3 at the same physical error rate.
        let p = 0.003;
        let rate = |d: u32| memory_experiment(d, d, p, 20_000, 11).unwrap().per_round();
        let (r3, r5) = (rate(3), rate(5));
        assert!(r5 < r3, "d=3 {r3}, d=5 {r5}");
    }
}