Swap γ and δ usage in the Clifford construction, Fix the Pauli‐gate mapping, and Use circ.add_circuit in‐place

[AiPals]

Swap γ and δ usage in the Clifford construction so that
γ → CX and S (diagonal)
δ → CZ
The original had these reversed.
Fix the Pauli‐gate mapping:
0 → Identity (no gate)
1 → X
2 → Y
3 → Z
Use circ.add_circuit in‐place rather than reassigning its (None) return value.
Cover all four Pauli values and skip the identity.

Code:

I have identify 4 errors on your code.

Fixes by Jesus Carrasco

from typing import List, Tuple

import numpy as np
from pytket import Circuit

def sample_q_mallows(n_qubits: int) -> Tuple[List[int], List[int]]:
hadamard_layer = [0] * n_qubits
permutation = [0] * n_qubits
avail_qubits = list(range(n_qubits))
log2 = np.log(2.0)

for i in range(n_qubits):
    m = len(avail_qubits)
    r = np.random.uniform(0, 1)
    index = int(2*m - np.ceil(np.log(r*(4**m - 1) + 1) / log2))
    hadamard_layer[i] = int(index < m)
    k = index if index < m else (2*m - index - 1)
    permutation[i] = avail_qubits[k]
    del avail_qubits[k]

return hadamard_layer, permutation

def clifford_canonical_F(
pauli_layer: List[int],
gamma: np.ndarray,
delta: np.ndarray
) -> Circuit:
"""Constructs an H-free Clifford in O·P·CZ·CX order (with γ→CX,S; δ→CZ)."""
n = len(pauli_layer)
circ = Circuit(n)

# 1) CX layer from gamma
for j in range(n):
    for i in range(j):
        if gamma[i, j]:
            circ.CX(i, j, opgroup="Clifford 2")

# 2) CZ layer from delta
for j in range(n):
    for i in range(j):
        if delta[i, j]:
            circ.CZ(i, j, opgroup="Clifford 2")

# 3) S-gates on delta diagonal
for i in range(n):
    if delta[i, i]:
        circ.S(i, opgroup="Clifford 1")

# 4) Pauli layer: 1→X, 2→Y, 3→Z, 0→I (skip)
for i, gate in enumerate(pauli_layer):
    if gate == 1:
        circ.X(i, opgroup="Clifford 1")
    elif gate == 2:
        circ.Y(i, opgroup="Clifford 1")
    elif gate == 3:
        circ.Z(i, opgroup="Clifford 1")
    # gate == 0: identity—no action

return circ

def find_random_gamma_delta(
n_qubits: int,
hadamard_layer: List[int],
permute_layer: List[int]
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""Generates two (γ, δ) pairs for H-free Cliffords in the canonical form."""
Delta1 = np.eye(n_qubits, dtype=int)
Delta2 = Delta1.copy()
Gamma1 = np.zeros((n_qubits, n_qubits), dtype=int)
Gamma2 = np.zeros_like(Gamma1)

# Diagonal bits
for i in range(n_qubits):
    Gamma2[i, i] = np.random.randint(2)
    if hadamard_layer[i]:
        Gamma1[i, i] = np.random.randint(2)

# Off‐diagonals
for i in range(n_qubits):
    for j in range(i + 1, n_qubits):
        b = np.random.randint(2)
        Gamma2[i, j] = Gamma2[j, i] = b
        Delta2[i, j] = np.random.randint(2)

        # Conditions for Gamma1
        if hadamard_layer[i] and hadamard_layer[j]:
            b = np.random.randint(2)
            Gamma1[i, j] = Gamma1[j, i] = b
        if hadamard_layer[i] and not hadamard_layer[j] and permute_layer[i] < permute_layer[j]:
            b = np.random.randint(2)
            Gamma1[i, j] = Gamma1[j, i] = b
        if not hadamard_layer[i] and hadamard_layer[j] and permute_layer[i] > permute_layer[j]:
            b = np.random.randint(2)
            Gamma1[i, j] = Gamma1[j, i] = b

        # Conditions for Delta1
        if not hadamard_layer[i] and hadamard_layer[j]:
            Delta1[i, j] = np.random.randint(2)
        if hadamard_layer[i] and hadamard_layer[j] and permute_layer[i] > permute_layer[j]:
            Delta1[i, j] = np.random.randint(2)
        if not hadamard_layer[i] and not hadamard_layer[j] and permute_layer[i] < permute_layer[j]:
            Delta1[i, j] = np.random.randint(2)

return Delta1, Delta2, Gamma1, Gamma2

def random_clifford_circ(n_qubits: int, **kwargs) -> Circuit:
"""
Samples a random n-qubit Clifford F·H·S·F' in canonical form using the quantum Mallows law.
"""
np.random.seed(kwargs.get("seed", None))
circ = Circuit(n_qubits)

hadamard, permute = sample_q_mallows(n_qubits)
D1, D2, G1, G2 = find_random_gamma_delta(n_qubits, hadamard, permute)

# First H-free Clifford F' (random Pauli)
pauli_rand = [np.random.randint(1, 4) for _ in range(n_qubits)]
Fp = clifford_canonical_F(pauli_rand, G2, D2)
circ.add_circuit(Fp, list(range(n_qubits)), [])

# Permutation layer S
for i in range(n_qubits):
    while permute[i] != i:
        j = permute[i]
        circ.SWAP(i, j, opgroup="Clifford 2")
        permute[i], permute[j] = permute[j], permute[i]

# Hadamard layer H
for i, h in enumerate(hadamard):
    if h:
        circ.H(i, opgroup="Clifford 1")

# Second H-free Clifford O (identity Pauli)
O = clifford_canonical_F([0]*n_qubits, G1, D1)
circ.add_circuit(O, list(range(n_qubits)), [])

return circ