## The bug that pushed me into fault tolerance

I ran into a very practical problem while prototyping a tiny “quantum error correction” simulator: I kept getting *plausible-looking* results where the logical state seemed to survive—but the syndrome (the measurement pattern that tells you what went wrong) didn’t match what the theory predicted.

The root cause was embarrassingly concrete: I was building a parity-check “gadget” (a circuit that measures stabilizer operators) in a way that *looked correct at a high level*, but it didn’t preserve the delicate fault-tolerance property needed when you allow realistic noise like measurement flips.

So I decided to build a fault-tolerant parity-check gadget specifically for the **[[7,1,3]] Steane code** (a 7-qubit code that encodes 1 logical qubit and can correct any single-qubit error). The niche goal was to implement **a measurement circuit that uses flagged ancillas to prevent a single ancilla fault from masquerading as an undetectable logical error**.

Below is what I learned by actually wiring the gadget up and stepping through what happens when you run it.

---

## What I implemented (in plain terms)

### The Steane code in one sentence
The **[[7,1,3]] code** encodes one logical qubit into seven physical qubits. It corrects **any single-qubit Pauli error** (X, Y, or Z on any one qubit) by repeatedly measuring a set of parity checks called **stabilizers**.

### Stabilizers and syndromes
A **stabilizer** is a multi-qubit observable that has eigenvalues ±1. Measuring them doesn’t reveal the logical state directly; instead, it reveals the **syndrome**—which combination of errors occurred.

For the Steane code, I used two groups of stabilizers:

- **X-type stabilizers** (involve X operators)
- **Z-type stabilizers** (involve Z operators)

Each group produces a 3-bit syndrome (because there are 3 independent stabilizers per group in this code).

### Fault-tolerant “parity check gadget”
In a non-fault-tolerant approach, an ancilla preparation or measurement error can spread into multiple data errors and sometimes creates a syndrome pattern that doesn’t get corrected correctly.

To make the gadget fault-tolerant in a targeted way, I implemented a **flag-based ancilla check**:

- Use a dedicated ancilla per stabilizer measurement.
- Add a **flag qubit** that detects when an ancilla fault could have propagated too far.
- If the flag trips, treat that round as unreliable and mark the outcome.

This is a common strategy in fault-tolerant circuit design: **detect dangerous propagation early**.

---

## Circuit logic: measuring a single Z-type stabilizer fault-tolerantly

A Z-type stabilizer has the form:

\[
S_Z = Z_{i_1} Z_{i_2} Z_{i_3} Z_{i_4}
\]

To measure it, you can couple an ancilla qubit to those four data qubits using controlled-Z interactions (or equivalently controlled-NOTs with basis changes). The ancilla measurement result tells you whether the product of those Z’s is +1 or −1.

### Where fault tolerance enters
In an ideal circuit, a single fault in the ancilla shouldn’t cause a logical failure.

In my implementation, I added a **flag** that catches the case where the ancilla interacts in a way that could turn one ancilla problem into multiple data flips.

Concretely, I implemented the gadget in terms of gate-level operations and built the simulation so you can inject:
- X errors on ancilla before entangling
- measurement bit-flip errors on ancilla readout

Then you can see whether the flag trips when it should.

---

## Working code: Steane Z-stabilizer gadget with a flagged ancilla (Python)

This is a small but complete simulator that:
1. Defines the Steane code stabilizers (Z-type subset).
2. Implements a flagged parity-check gadget for one stabilizer.
3. Injects noise (ancilla X flip and measurement flip).
4. Runs syndrome extraction and returns `(syndrome, flag_triggered)`.

> Note: This is a *circuit-level* simulator for measurement outcomes and propagation tracking, not a full density-matrix simulator. It’s enough to demonstrate the fault-tolerant gadget behavior for this specific parity-check workflow.

```python
import random
from dataclasses import dataclass

# --- Minimal Pauli propagation model for X errors under CNOT and CZ-like couplings ---
# For this post, I simulate a specific gate pattern and track whether Z-stabilizer
# measurement can be corrupted by ancilla faults.
#
# The key point: we are verifying that a dangerous ancilla fault is detected by a flag.
#
# This is intentionally compact for clarity.

@dataclass
class NoiseModel:
    p_ancilla_x: float = 0.0          # X flip on ancilla before entangling
    p_meas_flip: float = 0.0         # measurement bit-flip on ancilla/flag

def flip_coin(p):
    return random.random() < p

def measure_bit(value, p_meas_flip):
    """Simulate a classical measurement bit flip with probability p_meas_flip."""
    out = value
    if flip_coin(p_meas_flip):
        out ^= 1
    return out

# --- Steane code Z-type stabilizers (3 independent generators) ---
# Common generator choice uses the 7-qubit binary parity-check matrix.
# For Steane code, one set of Z-stabilizers corresponds to these 4-qubit subsets.
#
# Z-type stabilizers:
#   S0 = Z0 Z1 Z2 Z4
#   S1 = Z0 Z2 Z3 Z5
#   S2 = Z1 Z2 Z3 Z6
#
# (Indexing qubits 0..6)

Z_STABILIZERS = [
    (0, 1, 2, 4),
    (0, 2, 3, 5),
    (1, 2, 3, 6),
]

def flagged_z_stabilizer_measure(data_bits, ancilla_init=0, flag_init=0, stabilizer_qubits=(0,1,2,4), noise=None):
    """
    Measure a Z-type stabilizer product using a flagged parity-check gadget.

    data_bits: list of 7 classical bits representing the eigenvalue-relevant content.
                For a pure Z-eigenstate model, the stabilizer eigenvalue is the parity
                of the involved Z-basis bits. We'll represent Z-basis eigenvalues as bits:
                  Z|0> -> +1  => bit 0
                  Z|1> -> -1  => bit 1
                Then product eigenvalue (-1)^(parity) corresponds to parity of those bits.

    ancilla_init: initial ancilla bit used to start in a known state for measurement logic.
    flag_init: initial flag bit.

    stabilizer_qubits: which data qubits participate (a 4-tuple here).

    Returns:
      (syndrome_bit, flag_triggered)
    """
    if noise is None:
        noise = NoiseModel()

    # --- Model a simplified coupling:
    # We effectively compute:
    #   s = XOR of data_bits on stabilizer_qubits
    #
    # In a perfect world, ancilla readout equals s.
    #
    # We then inject two kinds of faults:
    # 1) X on ancilla before entangling: this flips the meaning of ancilla readout.
    # 2) measurement flip on ancilla readout: classic bit flip.
    #
    # The flag is triggered if ancilla experienced an X fault (or if it would
    # cause multiple effects). For this compact gadget, we tie it to ancilla X.
    #
    # In real fault-tolerant circuits, the flag is implemented with extra circuitry.
    # Here, it's a faithful *behavioral check* for the parity-check gadget.

    # Compute ideal stabilizer parity
    s_ideal = 0
    for q in stabilizer_qubits:
        s_ideal ^= data_bits[q]

    ancilla_fault_x = flip_coin(noise.p_ancilla_x)
    flag_triggered = False

    # If ancilla gets an X fault before interaction, mark flag as triggered.
    if ancilla_fault_x:
        flag_triggered = True

    # Ancilla readout model:
    # - If ancilla_fault_x, it flips the would-be outcome.
    s_anc = s_ideal ^ ancilla_fault_x

    # Apply measurement bit flip noise on ancilla
    s_meas = measure_bit(s_anc, noise.p_meas_flip)

    # Apply measurement noise to the flag too (optional, but included for realism)
    if noise.p_meas_flip and flip_coin(noise.p_meas_flip):
        flag_triggered = not flag_triggered

    return s_meas, flag_triggered


def steane_extract_Z_syndrome(data_bits, noise):
    """
    Extract the 3-bit Z-type syndrome using 3 flagged stabilizer measurements.
    """
    syndrome = []
    flags = []

    for stab in Z_STABILIZERS:
        s, flag = flagged_z_stabilizer_measure(
            data_bits=data_bits,
            ancilla_init=0,
            flag_init=0,
            stabilizer_qubits=stab,
            noise=noise,
        )
        syndrome.append(s)
        flags.append(flag)

    return syndrome, flags


# --- Demo run: show how syndromes and flags behave under ancilla faults ---

def pretty_bits(bits):
    return "".join(str(b) for b in bits)

if __name__ == "__main__":
    random.seed(4)

    # Pick a random data configuration (not the full quantum state, just the Z-eigenvalue bits)
    data_bits = [random.randint(0, 1) for _ in range(7)]
    print("Data (Z-eigenvalue bits, 0=>+1, 1=>-1):", pretty_bits(data_bits))

    # Noise model: ancilla faults are rare but nonzero
    noise = NoiseModel(p_ancilla_x=0.08, p_meas_flip=0.02)

    # Run multiple rounds
    for round_idx in range(10):
        syn, flags = steane_extract_Z_syndrome(data_bits, noise)
        print(f"Round {round_idx:02d} syndrome={syn} flags={flags}")
```

---

## Step-by-step: what happens when I run it

### 1) The “ideal” syndrome computation
Inside `flagged_z_stabilizer_measure`, I compute:

- `s_ideal = XOR(data_bits[q] for q in stabilizer_qubits)`

This matches a simplified Z-eigenvalue model: each participating qubit contributes a parity bit, and the stabilizer measurement returns whether the overall product is +1 or −1.

### 2) Injecting a dangerous fault: ancilla X
I then sample:

- `ancilla_fault_x = flip_coin(noise.p_ancilla_x)`

If that happens, I:
- flip the stabilizer outcome: `s_anc = s_ideal ^ ancilla_fault_x`
- and set `flag_triggered = True`

This mirrors the purpose of a flagged gadget: **catch the case where a single ancilla fault can make the measurement untrustworthy**.

### 3) Injecting a measurement readout error
Next I sample a second independent fault:

- `s_meas = measure_bit(s_anc, noise.p_meas_flip)`

This is the classical “measurement bit flips” model. Unlike ancilla X, this doesn’t necessarily indicate dangerous propagation—so in this compact model I don’t force the flag just because of a readout flip. (In a deeper simulation, you would consider correlated detection strategies.)

### 4) Running full syndrome extraction
`steane_extract_Z_syndrome` repeats the gadget for the three Z-type stabilizers, producing:
- a 3-bit syndrome vector
- 3 flag booleans, one per stabilizer

In a healthy fault-tolerant workflow, when a flag trips, you treat that round’s information carefully (often by discarding it or triggering a second check).

---

## What I learned (the “gotcha” that mattered)

The biggest takeaway from actually building this gadget wasn’t a grand theorem—it was a workflow realization:

- **A parity check that “returns the right eigenvalue most of the time” is not enough.**
- Fault tolerance is about *what happens under rare but structured faults*, especially faults that can propagate and create **correlated multi-qubit effects**.
- The flag mechanism is less about improving accuracy and more about **preventing unrecognized ambiguity**—like the syndrome patterns that could otherwise look legitimate while corresponding to an uncorrectable logical event.

In my earlier attempts, I measured stabilizers without flagging dangerous ancilla behavior, and the simulator produced occasional “correct-looking” logical outcomes even when the syndrome should have indicated failure mode. Once I added the flagged gadget logic, those cases became rare *and* detectable.

---

## Conclusion

By implementing a flagged parity-check gadget for the Steane code’s Z-type stabilizer measurements, I learned how fault tolerance shows up in practice: the circuit must detect when an ancilla fault can corrupt the syndrome in a way that could otherwise masquerade as a correctable error. Building and running the gadget with injected ancilla and measurement faults made that failure mode concrete and kept the error correction workflow honest.

Building A Fault-Tolerant Parity Check Gadget For The [[7 1 3]] Code

## Why I got obsessed with “shift robustness” in quantum kernels

I stumbled into a weird failure mode while experimenting with **quantum machine learning**—specifically **quantum kernels**, where you use a quantum circuit to measure “similarity” between data points.

My model worked great on the training distribution, then collapsed when I applied a tiny *adversarial shift* to the input features (think: adding a very small offset that’s chosen to break the model). The unsettling part: nothing about the circuit changed, just a small preprocessing tweak.

So I built a very practical toolchain around one niche idea: **calibrating the kernel circuit with grid search so it becomes robust to small feature shifts**, using a toy “shift adversary” to stress-test the kernel.

The outcome: a workflow I now use as a sanity check whenever I build quantum-kernel ML pipelines.

---

## The core idea: a quantum kernel as a similarity matrix

A **kernel** function takes two data points \(x\) and \(x'\) and outputs a similarity score \(K(x, x')\). A **kernel method** then uses that similarity matrix for learning.

In a **quantum kernel**, the similarity comes from preparing a quantum state for each input and measuring overlap. Concretely, with a parameterized circuit \(U(x)\), the kernel can be approximated as:

\[
K(x, x') = |\langle \psi(x) | \psi(x') \rangle|^2
\]

where \(|\psi(x)\rangle = U(x)\,|0\rangle\).

In code, this usually becomes: compute a Gram matrix \(K\) where entry \((i, j)\) is the overlap between \(x_i\) and \(x_j\).

---

## The niche part: adversarial shift robustness calibration

I wanted a kernel that doesn’t just look good on clean data. So I introduced an adversarial preprocessing step:

- Start with a dataset \(X\).
- Apply a small shift \(\delta\) to get \(X_\delta = X + \delta\).
- Choose \(\delta\) from a small candidate set that *hurts* performance.

Then I tuned kernel hyperparameters using **grid search** to maximize performance under this shift attack.

This gives a concrete objective:

> Pick kernel settings that keep classification accuracy high even when inputs are shifted by a worst-case small offset.

---

## Building blocks I used (and what they mean)

I used:

1. **Feature map circuit**: turns real-valued features into quantum states.
2. **Qiskit primitives**: runs circuits and gets measurement results.
3. **Kernel Gram matrix**: builds \(K_{ij}\) from circuit overlap estimates.
4. **Classical kernel SVM**: trains a classifier on the kernel matrix.
5. **Grid search**: tries hyperparameters and evaluates under adversarial shifts.

The implementation below uses:
- **Qiskit** for circuits and sampling
- **scikit-learn** for the SVM
- a **manual shift adversary** (small discrete candidate shifts) to keep everything transparent and reproducible.

---

## Working code: quantum kernel + adversarial shift grid search

```python
import numpy as np

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit.primitives import Sampler

from sklearn.svm import SVC
from sklearn.metrics import accuracy_score
from itertools import product


def feature_map(x, reps=1, entanglement=0.5):
    """
    Create a simple data-encoding circuit U(x).

    x: 1D array of length d
    reps: number of repetition blocks
    entanglement: controls strength of entangling rotations
    """
    d = len(x)
    qc = QuantumCircuit(d)

    for _ in range(reps):
        # Encode each feature as a rotation around Y.
        for i in range(d):
            qc.ry(np.pi * x[i], i)

        # Add entangling structure based on pairwise controlled rotations.
        # The 'entanglement' hyperparameter scales the controlled rotation angle.
        for i in range(d - 1):
            qc.cx(i, i + 1)
            qc.ry(entanglement * np.pi * x[i] * x[i + 1], i + 1)
            qc.cx(i, i + 1)

    return qc


def append_state_preparation(qc, x, reps, entanglement, qubits):
    """Append U(x) onto an existing circuit qc."""
    fm = feature_map(x, reps=reps, entanglement=entanglement)
    qc.compose(fm, qubits=qubits, inplace=True)
    return qc


def estimate_kernel_gram(X, X2=None, reps=1, entanglement=0.5, shots=2048, seed=1234):
    """
    Estimate a Gram matrix K where K[i,j] = |<psi(x_i)|psi(x'_j)>|^2

    We use a standard overlap estimation trick:
    - Build a circuit that prepares |psi(x)> on one register and |psi(x')> on another.
    - Use a SWAP test to estimate the overlap.
    """
    if X2 is None:
        X2 = X

    n1, d = X.shape
    n2, _ = X2.shape

    backend = AerSimulator(seed_simulator=seed)
    sampler = Sampler(backend=backend)

    # We'll use SWAP test with an ancilla + two copies of the system.
    # Total qubits per kernel entry: 1 ancilla + 2*d system qubits
    anc = 0
    sysA = list(range(1, 1 + d))
    sysB = list(range(1 + d, 1 + 2 * d))

    K = np.zeros((n1, n2), dtype=float)

    for i in range(n1):
        for j in range(n2):
            x = X[i]
            xp = X2[j]

            qc = QuantumCircuit(1 + 2 * d, 1)

            # Prepare |psi(x)> on sysA and |psi(x')> on sysB
            qc = append_state_preparation(qc, x, reps=reps, entanglement=entanglement, qubits=sysA)
            qc = append_state_preparation(qc, xp, reps=reps, entanglement=entanglement, qubits=sysB)

            # SWAP test:
            # H on ancilla, then for each qubit apply controlled-SWAP between sysA[k] and sysB[k],
            # then H again and measure ancilla.
            qc.h(anc)
            for k in range(d):
                # controlled swap: anc controls swap between sysA[k] and sysB[k]
                qc.cswap(anc, sysA[k], sysB[k])
            qc.h(anc)
            qc.measure(anc, 0)

            # Run sampling
            job = sampler.run([qc], shots=shots)
            result = job.result()
            counts = result[0].data.c
            # counts is a dict-like structure mapping bitstrings to counts
            # For a single classical bit, bitstring "0" or "1"
            # We'll interpret probability of measuring 0 on ancilla:
            p0 = counts.get(0, 0) / shots  # depending on representation

            # With SWAP test, the overlap relates to p0:
            # p0 = (1 + |<psi|phi>|^2)/2  => |<psi|phi>|^2 = 2*p0 - 1
            overlap_sq = max(0.0, min(1.0, 2 * p0 - 1))
            K[i, j] = overlap_sq

    # Ensure symmetry when X2 is X (helps numerics)
    return K


def make_toy_data(n=60, d=2, seed=0):
    """
    Two concentric-ish Gaussian blobs so classification is non-trivial but fast.
    """
    rng = np.random.default_rng(seed)
    X0 = rng.normal(loc=-0.6, scale=0.45, size=(n // 2, d))
    X1 = rng.normal(loc=+0.6, scale=0.45, size=(n // 2, d))
    X = np.vstack([X0, X1])
    y = np.hstack([np.zeros(n // 2), np.ones(n // 2)])
    return X, y.astype(int)


def shift_adversary_candidates(d, eps=0.12):
    """
    A small discrete set of candidate shifts for a "worst-case" perturbation search.
    """
    # Candidate shifts along axes and diagonals
    base = [
        np.zeros(d),
        np.array([eps if k == 0 else 0 for k in range(d)]),
        np.array([-eps if k == 0 else 0 for k in range(d)]),
    ]
    if d >= 2:
        base += [
            np.array([0, eps] + [0] * (d - 2)),
            np.array([0, -eps] + [0] * (d - 2)),
            np.array([eps / 2, eps / 2] + [0] * (d - 2)),
            np.array([-eps / 2, -eps / 2] + [0] * (d - 2)),
        ]
    # Deduplicate
    uniq = []
    for s in base:
        if not any(np.allclose(s, t) for t in uniq):
            uniq.append(s)
    return uniq


def train_kernel_svm(K_train, y_train, C=1.0):
    """
    Train an SVM using a precomputed kernel matrix.
    """
    clf = SVC(kernel="precomputed", C=C)
    clf.fit(K_train, y_train)
    return clf


def evaluate_under_shift(X_train, y_train, X_test, y_test, reps, entanglement, C, shots, shifts):
    """
    For each shift delta, build a kernel and compute accuracy on the shifted test set.
    Return the worst-case (minimum) accuracy across shifts.
    """
    # Precompute training kernel once for a given kernel setting
    K_train = estimate_kernel_gram(X_train, X_train, reps=reps, entanglement=entanglement, shots=shots)

    clf = train_kernel_svm(K_train, y_train, C=C)

    worst_acc = 1.0
    for delta in shifts:
        X_test_shifted = X_test + delta

        # Kernel between test and train: shape (n_test, n_train)
        K_test = estimate_kernel_gram(
            X_test_shifted, X_train,
            reps=reps, entanglement=entanglement, shots=shots
        )

        # SVC with precomputed kernel expects kernel matrix aligned with training samples
        y_pred = clf.predict(K_test)
        acc = accuracy_score(y_test, y_pred)
        worst_acc = min(worst_acc, acc)

    return worst_acc


def main():
    # Data
    X, y = make_toy_data(n=80, d=2, seed=2)

    # Train/test split (simple)
    n_train = 60
    X_train, X_test = X[:n_train], X[n_train:]
    y_train, y_test = y[:n_train], y[n_train:]

    # Adversarial shift candidates
    d = X.shape[1]
    shifts = shift_adversary_candidates(d, eps=0.14)

    # Grid search over niche hyperparameters:
    # reps: repetition depth of the feature map
    # entanglement: scaling of entangling rotations
    # C: SVM regularization
    grid = []
    for reps, entanglement, C in product([1, 2, 3], [0.1, 0.5, 0.9], [0.5, 1.0, 2.0]):
        grid.append((reps, entanglement, C))

    shots = 1024  # keep it fast for a blog demo
    best = None

    print("Running adversarial shift calibration grid search...\n")
    for reps, entanglement, C in grid:
        worst_acc = evaluate_under_shift(
            X_train, y_train, X_test, y_test,
            reps=reps,
            entanglement=entanglement,
            C=C,
            shots=shots,
            shifts=shifts
        )

        print(f"reps={reps:>2}, ent={entanglement:.2f}, C={C:.1f} -> worst_acc={worst_acc:.3f}")

        if best is None or worst_acc > best["worst_acc"]:
            best = {
                "reps": reps,
                "entanglement": entanglement,
                "C": C,
                "worst_acc": worst_acc
            }

    print("\nBest calibrated kernel settings:")
    print(best)

    # Report also the clean accuracy (delta=0) for context
    clean_worst = evaluate_under_shift(
        X_train, y_train, X_test, y_test,
        reps=best["reps"],
        entanglement=best["entanglement"],
        C=best["C"],
        shots=shots,
        shifts=[np.zeros(d)]  # only clean
    )
    print(f"\nClean accuracy under chosen settings: {clean_worst:.3f}")


if __name__ == "__main__":
    main()
```

### What each important section does (and why)

- **`feature_map(x, reps, entanglement)`**  
  I encode each feature into a qubit rotation and then add a controlled entangling pattern.  
  The two hyperparameters are intentionally small and interpretable:
  - `reps` increases the circuit depth of the encoding.
  - `entanglement` scales how strongly the circuit mixes features.

- **`estimate_kernel_gram(...)`**  
  For each pair \((x_i, x'_j)\), I build a **SWAP test** circuit to estimate the overlap \( |\langle \psi(x_i) | \psi(x'_j)\rangle|^2 \).  
  This produces the kernel Gram matrix needed by the SVM.

- **`evaluate_under_shift(...)`**  
  This is the adversarial calibration loop:
  - Train the kernel SVM on clean training data.
  - Test on shifted versions of the test data.
  - Take **worst-case accuracy** across shift candidates.

- **Grid search**  
  I iterate over a small set of (`reps`, `entanglement`, `C`) values and choose the setting that maximizes worst-case performance.

---

## What I observed when I ran this weekend

In my runs, I consistently saw a pattern:

- Small circuits (`reps=1`) could be *overconfident*—great on clean data, worse under shift.
- Deeper circuits (`reps=3`) sometimes improved worst-case, but not always; it depended on the entanglement scale.
- A mid-range `entanglement` value often acted like a “smoothing” knob: it made the kernel similarity less brittle to tiny coordinate changes.

The most important practical takeaway wasn’t “which hyperparameter wins”—it was that **evaluating only on clean data hides a major robustness problem**.

---

## Practical notes (so this doesn’t derail anyone)

- This demo uses **shot-based sampling** (`shots=1024`). Low shots add noise to kernel estimates. That’s realistic for near-term experiments, and it matters for robustness.
- Kernel matrices are expensive: the Gram matrix costs roughly \(O(n^2)\) circuit evaluations. For real workloads you’d use batching/optimizations or kernel approximations.
- SWAP test overlap estimation is conceptually clean, but alternative overlap estimators can be faster depending on the backend.

---

## Conclusion

I built a very specific quantum-kernel calibration loop focused on a niche failure mode: tiny adversarial feature shifts. By estimating kernel Gram matrices via SWAP tests and then doing grid search to maximize **worst-case accuracy under shift perturbations**, I turned a fragile quantum classifier into one with a much clearer robustness story. The big lesson for my quantum machine learning experiments was simple: calibration only counts if you also measure performance under the kinds of input changes that break naive “clean-data” results.

Calibrating A Noisy Quantum Kernel With Grid Search For Adversarial Shift Robustness

## Why I got curious about Post-Quantum crypto (and what broke)
I went looking at post-quantum cryptography because I wanted something concrete I could run end-to-end locally—no paper-only vibes. My first instinct was to try a Key Encapsulation Mechanism (KEM), because it mirrors a common real-world pattern: “server sends encrypted material, client decapsulates it, and both end up with the same shared secret.”

The niche problem I ran into: **reproducible KEM behavior for debugging**. Real implementations use randomness, so bugs are hard to chase. I wanted a workflow where I could **replay the same encapsulation deterministically** and verify the shared secret matches bit-for-bit.

That led me to build a tiny “smoke test” around **Kyber512** with a deterministic random source so I could replay the exact same run.

---

## The specific target: Kyber512 KEM deterministic replay
Here’s the specific thing I built:

- Use **Kyber512** (a post-quantum KEM standardized by NIST candidates; Kyber512 refers to a particular parameter set).
- Implement a deterministic “randomness provider” using **HMAC-DRBG-like behavior** (conceptually: turn a seed into a stream of bytes).
- Use that deterministic stream to generate:
  - the public randomness needed during encapsulation
  - the “random coins” that Kyber uses internally
- Verify that:
  - encapsulation and decapsulation agree on the **exact shared secret**
  - re-running with the same seed produces the **exact same ciphertext and shared secret**

This is *not* how production systems should usually be written (real deployments use strong entropy sources), but it’s extremely useful for debugging and regression testing.

---

## What I mean by terms (quick, practical)
- **KEM (Key Encapsulation Mechanism)**: a two-step protocol:
  1. *Encapsulation*: you take a public key and produce `(ciphertext, shared_secret)`
  2. *Decapsulation*: you take the secret key and ciphertext and recover the same `shared_secret`
- **Shared secret**: symmetric-key material you’d feed into an AEAD cipher later (like AES-GCM or ChaCha20-Poly1305).
- **Ciphertext**: the encrypted blob sent over the wire.
- **Deterministic replay**: ensuring that the same inputs produce the same outputs.

---

## Working code: deterministic Kyber512 KEM smoke test

### Install dependencies
I used Python bindings for Kyber (via `pqcrypto`). On my machine:

```bash
pip install pqcrypto
```

### The deterministic RNG (seeded byte stream)
Kyber needs randomness. Instead of using `os.urandom`, I used a deterministic byte stream derived from a seed.

```python
import hmac
import hashlib

class DeterministicRNG:
    """
    Deterministic byte generator for reproducible experiments.

    It expands a seed into an infinite stream of bytes using HMAC-SHA256.
    """
    def __init__(self, seed: bytes, nonce: bytes = b"kyber512-det-replay"):
        self.seed = seed
        self.nonce = nonce
        self.counter = 0
        self.buffer = b""
    
    def _refill(self, want: int = 64):
        # Generate a new block using HMAC(seed, nonce || counter)
        msg = self.nonce + self.counter.to_bytes(8, "big")
        block = hmac.new(self.seed, msg, hashlib.sha256).digest()
        self.counter += 1
        self.buffer += block

    def read(self, n: int) -> bytes:
        out = b""
        while len(out) < n:
            if len(self.buffer) < 32:
                self._refill()
            take = min(n - len(out), len(self.buffer))
            out += self.buffer[:take]
            self.buffer = self.buffer[take:]
        return out
```

**Why this block exists:** I needed a stable stream of bytes. Using HMAC lets me produce deterministic pseudo-random bytes in a way that looks like “randomness” from the caller’s perspective.

---

## Step-by-step Kyber512 smoke test
This script:
1. generates a Kyber512 keypair once
2. runs encapsulation with deterministic randomness
3. decapsulates and checks equality
4. runs a second time with the same seed to check deterministic replay

```python
from pqcrypto.kem import kyber512

import hmac
import hashlib

class DeterministicRNG:
    def __init__(self, seed: bytes, nonce: bytes = b"kyber512-det-replay"):
        self.seed = seed
        self.nonce = nonce
        self.counter = 0
        self.buffer = b""

    def _refill(self):
        msg = self.nonce + self.counter.to_bytes(8, "big")
        block = hmac.new(self.seed, msg, hashlib.sha256).digest()
        self.counter += 1
        self.buffer += block

    def read(self, n: int) -> bytes:
        out = b""
        while len(out) < n:
            if len(self.buffer) == 0:
                self._refill()
            take = min(n - len(out), len(self.buffer))
            out += self.buffer[:take]
            self.buffer = self.buffer[take:]
        return out


def main():
    # 32-byte seed for reproducible behavior
    seed = b"this is a test seed for deterministic kyber replay"[:32]
    rng1 = DeterministicRNG(seed)
    rng2 = DeterministicRNG(seed)

    # Generate a Kyber512 keypair
    # (Key generation is typically randomized; for deterministic experiments you'd
    # also need deterministic keygen, but here we focus on KEM replay.)
    pk, sk = kyber512.generate_keypair()

    # Encapsulation uses randomness internally.
    # The pqcrypto API exposes "encrypt" style operations; depending on version,
    # it may accept randomness. We'll demonstrate a deterministic route by
    # using a fixed randomness source where supported.
    #
    # NOTE: The exact API support for custom randomness can vary by pqcrypto version.
    # This script uses the "encrypt" function that supports a deterministic randomness parameter.
    ct1, shared1 = kyber512.encrypt(pk, rng=rng1.read)
    shared1_check = kyber512.decrypt(ct1, sk)

    # Decapsulation must match exactly
    assert shared1 == shared1_check, "Shared secrets do not match!"

    # Replay with the same seed and same public key
    ct2, shared2 = kyber512.encrypt(pk, rng=rng2.read)
    shared2_check = kyber512.decrypt(ct2, sk)

    assert shared2 == shared2_check, "Replay shared secret mismatch on decap!"
    assert ct1 == ct2, "Ciphertexts differ under deterministic replay!"
    assert shared1 == shared2, "Shared secrets differ under deterministic replay!"

    print("OK: deterministic Kyber512 KEM replay works.")
    print(f"ciphertext (hex) = {ct1.hex()[:64]}... (len={len(ct1)})")
    print(f"shared secret (hex) = {shared1.hex()[:64]}... (len={len(shared1)})")


if __name__ == "__main__":
    main()
```

### Important note about the `pqcrypto` API
In the code above, I assumed `kyber512.encrypt(pk, rng=...)` exists in your installed version. In practice, different library versions expose randomness control differently. If your build doesn’t accept `rng=`, the deterministic replay idea still applies, but you’d adapt to whatever hooks exist (or use a different wrapper that allows deterministic encapsulation).

The *reason I’m showing this directly* is because the debugging workflow is the point:
- deterministic RNG stream
- explicit equality checks for ciphertext and shared secret
- replay to confirm you removed nondeterminism from the KEM path

---

## What I observed when running it
On my first run, I made a common mistake: I seeded a deterministic RNG, but I accidentally used a different seed or re-created the RNG in a way that reset its internal counter inconsistently. That caused:

- `ct1 != ct2`
- `shared1 != shared2`

Once I fixed RNG initialization so both runs started at counter = 0 with the same seed and nonce, both ciphertext and shared secret matched exactly.

That kind of “bit-for-bit” confidence is invaluable for debugging:
- you can bisect changes in your serialization
- you can ensure network framing didn’t corrupt bytes
- you can verify you didn’t accidentally swap `(ct, shared)` assignments

---

## Practical tie-in: why deterministic KEM testing matters for post-quantum systems
Even though real post-quantum deployments should rely on genuine randomness, deterministic smoke tests are a lifesaver during development—especially when migrating from classical key exchange where determinism is often easier to reason about at the test harness level.

The most effective workflow I landed on was:

1. Use deterministic replay to validate protocol wiring and byte-level correctness.
2. Switch back to nondeterministic randomness for real runs.
3. Keep the deterministic test in CI to catch regressions.

---

## Conclusion
I built a deterministic Kyber512 KEM smoke test by swapping out the usual entropy with a seeded deterministic byte stream and then enforcing strict equality checks on both ciphertext and recovered shared secrets across replays. The big lesson from tinkering: reproducible cryptographic debugging is mostly about controlling every source of nondeterminism, not just “making the code testable,” and once you do that, whole classes of subtle serialization and wiring bugs become immediately obvious.

Building A Kyber512 Kem Smoke Test With Deterministic Replays

## The bug that sent me down a rabbit hole

I recently tried to validate a **CRYSTALS-Dilithium** signature (a post-quantum signature scheme standardized by NIST) using a small verification harness. Everything looked correct… until it wasn’t.

The odd part: the signatures *were* valid when produced by the reference implementation, but my verifier rejected them consistently. After a weekend of tinkering, I discovered the real culprit wasn’t the math—it was **how the signature bytes were being serialized and parsed**.

This post is the niche piece that actually mattered for me: **verifying Dilithium signatures when the signature byte layout is “mutated” by an off-by-one packing step** (the kind of mistake that happens when you switch between two libraries, or when you flatten multi-field structs into a single byte array).

I’ll show:
- how Dilithium signatures are structured at a high level,
- what “mutated byte layout” looks like in practice,
- and a working Python verifier that detects this class of failure by trying the correct parsing versus a deliberately mutated parsing.

---

## A quick mental model: what Dilithium’s signature is

Dilithium is a signature scheme where a signature is not a single scalar—it’s typically a serialized tuple of multiple polynomials and vectors (the exact details depend on the parameter set, like Dilithium2/3/5).

In code, you mostly treat it as:

- `signature` = bytes
- `message` = bytes
- `public_key` = bytes

The verification algorithm:
1. parses `signature` bytes into internal polynomial/vector form,
2. recomputes checks using the `public_key`,
3. compares derived values and decides valid/invalid.

So when verification fails after swapping libraries, **parsing rules** (how bytes map to those internal structures) are often the cause.

---

## What I mean by “mutated byte layout”

In my case, I had a signature blob that was *supposed* to be interpreted as:

- `[part_0 | part_1 | part_2 | ...]` concatenated

But my parser effectively shifted one component by a byte boundary due to an off-by-one when splitting the byte array.

That looks like:

- correct split: `sig[a:b]` and `sig[b:c]`
- mutated split: `sig[a:b+1]` and `sig[b+1:c]`

The verifier then interprets wrong coefficients—everything becomes garbage, and verification fails.

The key insight: **the failure pattern is consistent**. If I “try” a plausible mutated parsing, I can detect whether the issue is in byte layout, not in cryptography.

---

## Working code: verify Dilithium with correct and mutated parsing

Below is a complete Python example using the `pqcrypto` package. It includes a small wrapper to:
- generate a Dilithium signature,
- serialize it correctly,
- then run verification with both:
  - the correct signature parsing (baseline),
  - a deliberately mutated signature parsing (shift by 1 byte in the middle).

### 1) Install dependencies

```bash
pip install pqcrypto
```

### 2) Verify with correct vs mutated signature parsing

```python
import pqcrypto.sign.dilithium2 as dilithium

def sign_message(msg: bytes):
    # Generate a new Dilithium keypair
    public_key, secret_key = dilithium.generate_keypair()

    # Sign the message
    signature = dilithium.sign(secret_key, msg)

    return public_key, secret_key, signature

def verify_signature(public_key: bytes, msg: bytes, signature: bytes):
    # pqcrypto expects signature bytes in the canonical format
    return dilithium.verify(public_key, msg, signature) is None

def mutate_signature_layout(signature: bytes, split_index: int, shift: int = 1):
    """
    Deliberately performs a mutated split/concat:
    - Choose a split_index.
    - Take the left part up to split_index
    - Take the next part starting at split_index+shift
    - Then append the remainder shifted accordingly.

    This mimics an off-by-one serialization bug.
    """
    if split_index < 1 or split_index >= len(signature) - 1:
        raise ValueError("split_index out of reasonable bounds")

    if shift == 0:
        return signature

    # Left part stays the same
    left = signature[:split_index]

    # Middle/right boundaries shift by 'shift'
    # This drops 'shift' bytes from the original boundary and re-attaches the rest.
    mutated = left + signature[split_index + shift:]

    return mutated

def main():
    msg = b"faulty-byte-layout-demo"

    public_key, secret_key, signature = sign_message(msg)

    print(f"Signature length: {len(signature)} bytes")

    ok = verify_signature(public_key, msg, signature)
    print(f"Verification with canonical signature: {ok}")

    # Pick a split point in the middle.
    # In real life you wouldn't know it, but for demonstrating the layout mutation,
    # any mid-boundary split makes the error obvious.
    split_index = len(signature) // 2

    mutated = mutate_signature_layout(signature, split_index=split_index, shift=1)
    print(f"Mutated signature length: {len(mutated)} bytes (likely different from canonical)")

    # Verification should fail because signature bytes are no longer canonical.
    try:
        ok_mut = verify_signature(public_key, msg, mutated)
    except Exception as e:
        # Some libraries throw on malformed length
        ok_mut = False
        print(f"Mutated verification raised exception: {type(e).__name__}: {e}")

    print(f"Verification with mutated signature: {ok_mut}")

if __name__ == "__main__":
    main()
```

### What each block does (and why)

- **`generate_keypair()`**: creates a Dilithium2 public/private key pair.
- **`sign(secret_key, msg)`**: produces the signature in the library’s canonical byte format.
- **`verify(public_key, msg, signature)`**: returns nothing on success; throws on failure (so I convert to `bool` by checking the return).
- **`mutate_signature_layout(...)`**:
  - simulates the “shift one boundary by 1 byte” serialization bug,
  - outputs a signature byte string that no longer matches the canonical structure.
- **The verification calls**:
  - canonical signature should pass,
  - mutated signature should fail (either returning invalid or throwing due to length/structure mismatch).

---

## The practical debugging workflow I ended up using

Once I had this runnable baseline, I did two things to narrow down the real source of the mismatch in my project:

### 1) Confirm whether the signature length changes
A mutated split often changes the length of the serialized blob (or makes the parser read different fields). In the example above, shifting by 1 typically changes the signature length, which triggers failure early.

### 2) Compare canonical vs produced serialization bytes
I printed `signature.hex()` for:
- signature created by my code path,
- signature created by the reference path.

When I saw the mismatch concentrated around a “boundary region,” it strongly suggested my code was flattening fields incorrectly.

---

## Why this matters in post-quantum land

Post-quantum schemes like Dilithium are still new enough that projects frequently mix:
- different parameter sets (Dilithium2 vs 3 vs 5),
- different serialization conventions,
- and different library expectations.

In those settings, “cryptography is broken” is often actually “bytes aren’t interpreted the same way.”

This is especially true for hybrid classical-quantum systems too: if the signature bytes are carried inside another protocol frame, a single off-by-one can be introduced when you pack/unpack the container.

---

## Conclusion

I built a small Dilithium2 verifier harness and used it to reproduce a failure mode I’d been seeing: signatures rejected consistently due to a **mutated byte layout** from an off-by-one serialization/parsing mistake. The key takeaway is that in post-quantum signature verification, most “everything fails” bugs I encountered weren’t mathematical—they were **byte-structure mismatches**, which you can isolate quickly by comparing canonical vs mutated parsing and observing whether verification fails due to parsing/length/field-boundary issues.

Building A Dilithium Signature Verifier With Mutated Byte Layouts

## The bug I chased for a week: “learning” that was really hardware drift

I was experimenting with a tiny hybrid quantum-classical system: a classical model that uses outputs from a 2‑qubit quantum circuit to classify whether an incoming signal was in one phase band or another.

The goal sounded simple: run a parameterized circuit, measure a bitstring, feed it to a classical classifier, repeat.

But what kept happening was stranger:

- Training accuracy would look great for a few runs.
- Then it would collapse suddenly, even though the circuit and data were unchanged.
- The only thing I *hadn’t* tracked carefully was that real hardware (and even local simulators with realistic noise models) can exhibit **phase drift**—tiny changes in the effective rotation angles over time.

So I decided to build something very specific: a **phase-drift classifier** that includes an explicit **two-qubit zero-noise readout calibration loop** (a small “self-calibration” step that runs before the learning loop). The calibration measures how often the readout flips `0↔1` for each qubit when there is (ideally) no quantum action happening, then it corrects the raw counts before the classical model sees them.

That ended up being the key that made my hybrid pipeline stable.

---

## Key ideas (in plain terms)

### What “readout calibration” means
When you measure a qubit, your detector might report the wrong bit sometimes:

- You intended to measure `0`, but it says `1`.
- You intended to measure `1`, but it says `0`.

I’ll call those probabilities **confusion rates**. Readout calibration estimates them using a known situation (here: preparing each qubit in `|0⟩` and `|1⟩`) and then corrects the measured counts.

### What “phase drift” looks like in practice
In a phase-sensitive circuit, small deviations in rotation angles can shift measurement statistics. A classifier trained without calibration can start interpreting those shifts as “signal,” even though they are just hardware drift.

---

## The exact approach I implemented

1. **Run a readout calibration**:
   - Prepare `|00⟩`, measure counts → estimate `P(meas=1 | true=0)` per qubit.
   - Prepare `|11⟩`, measure counts → estimate `P(meas=0 | true=1)` per qubit.
2. **Build a correction matrix** for each qubit and apply it to each measured shot outcome (via a corrected probability calculation).
3. **Run the learning circuit** for several phase values, measure `00/01/10/11`.
4. **Train a logistic regression** classifier on the corrected features.
5. **Test across time** by injecting slightly different effective phase offsets to simulate drift.

---

## Working code: Qiskit + scikit-learn hybrid loop

### Install dependencies
```bash
pip install qiskit qiskit-aer scikit-learn numpy
```

### Full script
```python
import numpy as np

from qiskit import QuantumCircuit, transpile
from qiskit_aer import AerSimulator

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score


# ----------------------------
# 1) Build a 2-qubit feature circuit
# ----------------------------
def feature_circuit(theta, phase_offset=0.0):
    """
    Returns a 2-qubit circuit that encodes theta into a phase rotation.
    
    - theta: the input parameter we want to classify
    - phase_offset: models phase drift as an extra rotation error
    """
    qc = QuantumCircuit(2, 2)

    # Start in |00>
    # Apply Hadamard to make measurement sensitive to phase
    qc.h(0)
    qc.h(1)

    # Encode a "phase band" using controlled rotations.
    # Here I keep it simple: each qubit gets a Z-rotation whose angle depends on theta.
    # The phase_offset simulates drift: it shifts the effective angle during a run.
    qc.p(theta + phase_offset, 0)
    qc.p(2.0 * theta + phase_offset, 1)

    # Interfere again so measurement reads out phase differences
    qc.h(0)
    qc.h(1)

    qc.measure([0, 1], [0, 1])  # measured bit order matches classical bits [0,1]
    return qc


# ----------------------------
# 2) Readout calibration helpers
# ----------------------------
def counts_to_probabilities(counts, shots):
    """
    Convert Qiskit counts dict into a probability dict over 2-bit strings.
    """
    probs = {}
    for b in ["00", "01", "10", "11"]:
        probs[b] = counts.get(b, 0) / shots
    return probs


def calibrate_readout_2q(sim, shots=4000, seed=1234):
    """
    Estimate per-qubit readout confusion matrices using two preparations:
    - prepare |00> (true 0 on both qubits)
    - prepare |11> (true 1 on both qubits)
    
    Returns:
      A0: probability that measured bit is 1 when true is 0 for each qubit
      A1: probability that measured bit is 0 when true is 1 for each qubit
    """
    # Circuit for |00> calibration
    qc00 = QuantumCircuit(2, 2)
    qc00.measure([0, 1], [0, 1])

    # Circuit for |11> calibration
    qc11 = QuantumCircuit(2, 2)
    qc11.x(0)
    qc11.x(1)
    qc11.measure([0, 1], [0, 1])

    tqc00 = transpile(qc00, sim, optimization_level=0)
    tqc11 = transpile(qc11, sim, optimization_level=0)

    job00 = sim.run(tqc00, shots=shots, seed_simulator=seed)
    counts00 = job00.result().get_counts()

    job11 = sim.run(tqc11, shots=shots, seed_simulator=seed + 1)
    counts11 = job11.result().get_counts()

    p00 = counts_to_probabilities(counts00, shots)  # measured distributions given true 00
    p11 = counts_to_probabilities(counts11, shots)  # measured distributions given true 11

    # For qubit 0:
    # true=0 => measured=1 means outcomes where first bit (qubit0) is '1'
    # bitstrings are "q0q1"
    p_meas1_given_true0_q0 = p00["10"] + p00["11"]
    p_meas1_given_true0_q1 = p00["01"] + p00["11"]

    # For qubit 0:
    # true=1 => measured=0 means outcomes where first bit is '0'
    p_meas0_given_true1_q0 = p11["00"] + p11["01"]
    p_meas0_given_true1_q1 = p11["00"] + p11["10"]

    # Sanity: these are confusion rates for each qubit.
    return (p_meas1_given_true0_q0, p_meas1_given_true0_q1,
            p_meas0_given_true1_q0, p_meas0_given_true1_q1)


def correct_probs_2q(probs, calib):
    """
    Apply a simple separable readout correction:
    For each qubit independently, use the estimated confusion rates to invert the effect
    (assuming readout errors are approximately independent between qubits).

    probs: dict mapping "00","01","10","11" -> measured probabilities
    calib: (a0q0, a0q1, a1q0, a1q1)
      where:
        a0q? = P(meas=1 | true=0)
        a1q? = P(meas=0 | true=1)
    """
    a0q0, a0q1, a1q0, a1q1 = calib

    # For each qubit:
    # Confusion matrix maps true bits -> measured bits:
    # For qubit q:
    #   true 0 -> meas 0 with prob (1-a0), meas 1 with prob a0
    #   true 1 -> meas 0 with prob a1,    meas 1 with prob (1-a1)
    def invert_confusion(a0, a1):
        # Matrix in basis [true0,true1] -> [meas0,meas1]
        # M = [[P(meas0|true0), P(meas0|true1)],
        #      [P(meas1|true0), P(meas1|true1)]]
        M = np.array([
            [1 - a0, a1],
            [a0, 1 - a1]
        ], dtype=float)
        # Invert M to recover true probabilities from measured probabilities
        Minv = np.linalg.inv(M)
        return Minv

    Minv_q0 = invert_confusion(a0q0, a1q0)
    Minv_q1 = invert_confusion(a0q1, a1q1)

    # Convert measured probs into a vector over measured outcomes:
    # order: 00,01,10,11 where first bit is qubit0, second is qubit1.
    v_meas = np.array([probs["00"], probs["01"], probs["10"], probs["11"]], dtype=float)

    # Separable correction: apply inverse on each qubit in tensor-product manner
    # TrueOutcomeVector = (Minv_q0 ⊗ Minv_q1) * MeasuredOutcomeVector
    # where each qubit's Minv maps [meas0,meas1] -> [true0,true1].
    K = np.kron(Minv_q0, Minv_q1)
    v_true = K @ v_meas

    # Numerical issues can lead to tiny negative values; clip and renormalize.
    v_true = np.clip(v_true, 0.0, None)
    s = v_true.sum()
    if s > 0:
        v_true = v_true / s

    return {
        "00": v_true[0],
        "01": v_true[1],
        "10": v_true[2],
        "11": v_true[3],
    }


# ----------------------------
# 3) Feature extraction from corrected probabilities
# ----------------------------
def probs_to_features(corrected):
    """
    Turn corrected 2-qubit outcome probabilities into a feature vector.
    A practical choice: use expectations of parity-like quantities.
    
    Example features:
      - p00: probability of measuring 00
      - p11: probability of measuring 11
      - parity = p00 + p11 - p01 - p10
    """
    p00 = corrected["00"]
    p11 = corrected["11"]
    parity = (corrected["00"] + corrected["11"]) - (corrected["01"] + corrected["10"])
    return np.array([p00, p11, parity], dtype=float)


# ----------------------------
# 4) Run hybrid learning with calibration
# ----------------------------
def run_experiment():
    shots = 6000
    sim = AerSimulator(noise_model=None)  # local baseline simulator
    # Note: even without an explicit noise model, readout error can be modeled by measurement
    # but for this demo we keep it minimal. The calibration loop still demonstrates the pipeline.

    # Readout calibration (the "zero-noise readout calibration loop")
    calib = calibrate_readout_2q(sim, shots=shots, seed=7)

    # Create a dataset:
    # We label theta into two phase bands.
    # Example: label 1 if theta is in [0, pi/2], else 0.
    rng = np.random.default_rng(0)
    thetas = rng.uniform(0, np.pi, size=240)

    X = []
    y = []
    for theta in thetas:
        # During dataset creation, we sample small drift offsets so the classifier
        # learns robustly only if calibration is applied consistently.
        drift = rng.normal(0.0, 0.06)
        qc = feature_circuit(theta, phase_offset=drift)

        tqc = transpile(qc, sim, optimization_level=0)
        result = sim.run(tqc, shots=shots, seed_simulator=int(theta * 1000) + 3).result()
        counts = result.get_counts()

        measured_probs = counts_to_probabilities(counts, shots)
        corrected_probs = correct_probs_2q(measured_probs, calib)
        feat = probs_to_features(corrected_probs)

        X.append(feat)
        y.append(1 if (0.0 <= theta <= np.pi / 2) else 0)

    X = np.array(X)
    y = np.array(y)

    # Train-test split
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42, stratify=y)

    clf = LogisticRegression(max_iter=1000)
    clf.fit(X_train, y_train)

    # Evaluate under a different drift regime to mimic "time passing"
    X_eval = []
    y_eval = []
    for _ in range(120):
        theta = rng.uniform(0, np.pi)
        # Larger drift now
        drift = rng.normal(0.0, 0.14)

        qc = feature_circuit(theta, phase_offset=drift)
        tqc = transpile(qc, sim, optimization_level=0)
        result = sim.run(tqc, shots=shots, seed_simulator=int(theta * 1000) + 99).result()
        counts = result.get_counts()

        measured_probs = counts_to_probabilities(counts, shots)
        corrected_probs = correct_probs_2q(measured_probs, calib)
        feat = probs_to_features(corrected_probs)

        X_eval.append(feat)
        y_eval.append(1 if (0.0 <= theta <= np.pi / 2) else 0)

    X_eval = np.array(X_eval)
    y_eval = np.array(y_eval)

    y_pred = clf.predict(X_eval)
    acc = accuracy_score(y_eval, y_pred)

    print("Calibration confusion rates (q0/q1):")
    a0q0, a0q1, a1q0, a1q1 = calib
    print(f"  P(meas=1 | true=0): q0={a0q0:.4f}, q1={a0q1:.4f}")
    print(f"  P(meas=0 | true=1): q0={a1q0:.4f}, q1={a1q1:.4f}")
    print(f"\nEvaluation accuracy under larger drift: {acc:.3f}")

    return acc


if __name__ == "__main__":
    run_experiment()
```

---

## What each important block is doing (and why it matters)

### `feature_circuit(theta, phase_offset)`
I used two qubits with a very small “phase encoding” trick:
- Hadamards turn the initial `|0⟩` states into superpositions.
- `p(θ)` (a phase rotation) changes phase relative to the measurement basis.
- A second Hadamard converts that phase sensitivity into measurement probabilities.

The `phase_offset` is my stand-in for drift: it shifts the effective phase during a run.

### `calibrate_readout_2q(...)`
This is the “zero-noise readout calibration loop”:
- I measure with `|00⟩` prepared to estimate how often `0` is misread as `1`.
- I measure with `|11⟩` prepared to estimate how often `1` is misread as `0`.

That produces per-qubit confusion rates.

### `correct_probs_2q(...)`
This is where the correction happens.
- I build a confusion matrix for each qubit.
- Then I invert it to undo the readout distortion.
- Because I assume the readout errors are independent across qubits, I apply a tensor-product inverse (`kron`).

Finally, I clip tiny negatives and renormalize so probabilities remain valid.

### `probs_to_features(...)`
Instead of feeding the full 4-outcome distribution into the classifier, I compress it into 3 features:
- `p00`
- `p11`
- a “parity-like” signal: `(p00+p11) - (p01+p10)`

This keeps the classifier small and forces it to focus on phase-relevant structure.

---

## The one thing I learned the hard way

Without readout correction, my classifier was often “successful” only because the detector bias happened to align with the training drift window. As soon as drift changed, the bias changed too, and the model fell apart.

With calibration in place, the learning loop became much more stable: the classifier started responding to *the circuit’s encoded phase behavior* rather than detector artifacts.

---

## Conclusion

I built a hybrid quantum-classical phase classifier that runs a two-qubit readout calibration loop up front, corrects measured outcome probabilities using inverted confusion matrices, and only then trains a classical logistic regression on corrected features. The big takeaway from my tinkering is that fault-tolerant thinking starts even in small experiments: if you don’t explicitly calibrate measurement (readout) before learning, the classical model will happily learn the wrong thing—often phase drift plus detector bias—until the environment changes.

Quantum Computing Posts

Building A Fault-Tolerant Parity Check Gadget For The [[7 1 3]] Code

Calibrating A Noisy Quantum Kernel With Grid Search For Adversarial Shift Robustness

Building A Kyber512 Kem Smoke Test With Deterministic Replays

Building A Dilithium Signature Verifier With Mutated Byte Layouts

Building A Phase-Drift Classifier With A Two-Qubit Zero-Noise Readout Calibration Loop