Skip to content

Quantum

Quantum-classical hybrid spiking layers: hardware backend bridge, hybrid quantum-spiking circuits, quantum error correction, noise modelling, and parameter-shift gradient optimization.

Hardware Bridge

sc_neurocore.quantum.hardware_bridge

Hardware Bridge for Quantum-Classical Hybrid execution.

This module provides the interface to offload the simulated quantum stochastic logic to actual quantum hardware via Qiskit, or to high-fidelity tensor-network simulators via PennyLane.

Usage::

from sc_neurocore.quantum.hardware_bridge import QuantumHardwareLayer
layer = QuantumHardwareLayer(n_qubits=4, backend_type="aer_simulator")
out_bits = layer.forward(input_bitstreams)

QuantumHardwareLayer dataclass

Executes a Quantum-Classical Hybrid Layer on Qiskit/PennyLane. Maps bitstream probability -> Qubit Rotation -> True Measurement.

Source code in src/sc_neurocore/quantum/hardware_bridge.py 
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
@dataclass
class QuantumHardwareLayer:
    """
    Executes a Quantum-Classical Hybrid Layer on Qiskit/PennyLane.
    Maps bitstream probability -> Qubit Rotation -> True Measurement.
    """

    n_qubits: int
    length: int = 1024
    backend_type: str = "aer_simulator"  # "aer_simulator", "pennylane.default.qubit", etc.
    _qiskit_simulator: Any = None
    _pennylane_dev: Any = None

    def __post_init__(self) -> None:
        if self.backend_type.startswith("aer") and not HAS_QISKIT:
            from sc_neurocore.exceptions import SCDependencyError

            raise SCDependencyError("Qiskit is required for aer_simulator backend.")
        if self.backend_type.startswith("pennylane") and not HAS_PENNYLANE:
            from sc_neurocore.exceptions import SCDependencyError

            raise SCDependencyError("PennyLane is required for pennylane backend.")

        if self.backend_type == "aer_simulator":
            self._qiskit_simulator = AerSimulator()
        elif self.backend_type.startswith("pennylane"):
            self._pennylane_dev = qml.device(
                "default.qubit", wires=self.n_qubits, shots=self.length
            )

    def forward(self, input_bitstreams: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
        """
        input_bitstreams: (n_qubits, length)
        Returns:
        output_bitstreams: (n_qubits, length)
        """
        p_in = np.mean(input_bitstreams, axis=1)
        theta = p_in * np.pi

        if self.backend_type == "aer_simulator":
            return self._run_qiskit(theta)
        elif self.backend_type.startswith("pennylane"):
            return self._run_pennylane(theta)
        else:
            raise ValueError(f"Unknown backend: {self.backend_type}")

    def _run_qiskit(self, theta: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
        """Runs the circuit on Qiskit AerSimulator for `length` shots."""
        qc = QuantumCircuit(self.n_qubits, self.n_qubits)

        # Apply Ry rotations based on theta
        for i in range(self.n_qubits):
            qc.ry(theta[i], i)

        qc.measure(range(self.n_qubits), range(self.n_qubits))

        # Run circuit for self.length shots
        compiled_circuit = transpile(qc, self._qiskit_simulator)
        job = self._qiskit_simulator.run(compiled_circuit, shots=self.length)
        result = job.result()
        counts = result.get_counts(compiled_circuit)

        # Reconstruct bitstreams from shot counts
        out_bits = np.zeros((self.n_qubits, self.length), dtype=np.uint8)
        current_idx = 0
        for bitstring, count in counts.items():
            # bitstring is like '0101' where index 0 is the last character in string
            for i in range(count):
                if current_idx < self.length:
                    for qubit_idx in range(self.n_qubits):
                        # Qiskit orders bitstrings right-to-left
                        bit_val = int(bitstring[self.n_qubits - 1 - qubit_idx])
                        # Invert because measurement logic expects |0> as 1
                        out_bits[qubit_idx, current_idx] = 1 - bit_val
                    current_idx += 1

        return out_bits

    def _run_pennylane(self, theta: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
        """Runs the circuit on PennyLane for `length` shots."""

        @qml.qnode(self._pennylane_dev)
        def circuit(angles: np.ndarray[Any, Any]) -> Any:
            for i in range(self.n_qubits):
                qml.RY(angles[i], wires=i)
            return qml.sample()

        # Returns shape: (shots, n_qubits)
        samples = circuit(theta)
        # Transpose to (n_qubits, shots) and invert so |0> -> 1
        res: np.ndarray[Any, Any] = (1 - samples).T.astype(np.uint8)
        return res

forward(input_bitstreams)

input_bitstreams: (n_qubits, length) Returns: output_bitstreams: (n_qubits, length)

Source code in src/sc_neurocore/quantum/hardware_bridge.py 
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
def forward(self, input_bitstreams: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
    """
    input_bitstreams: (n_qubits, length)
    Returns:
    output_bitstreams: (n_qubits, length)
    """
    p_in = np.mean(input_bitstreams, axis=1)
    theta = p_in * np.pi

    if self.backend_type == "aer_simulator":
        return self._run_qiskit(theta)
    elif self.backend_type.startswith("pennylane"):
        return self._run_pennylane(theta)
    else:
        raise ValueError(f"Unknown backend: {self.backend_type}")

Hybrid Layer

sc_neurocore.quantum.hybrid

QuantumStochasticLayer dataclass

Simulates a Quantum-Classical Hybrid Layer. Input bitstream probability -> Qubit Rotation -> Measurement Probability.

Mapping: p_in -> theta = p_in * pi P_out = |<0|Ry(theta)|0>|^2 = cos^2(theta/2) This non-linearity is useful for classification.

Source code in src/sc_neurocore/quantum/hybrid.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
@dataclass
class QuantumStochasticLayer:
    """
    Simulates a Quantum-Classical Hybrid Layer.
    Input bitstream probability -> Qubit Rotation -> Measurement Probability.

    Mapping: p_in -> theta = p_in * pi
    P_out = |<0|Ry(theta)|0>|^2 = cos^2(theta/2)
    This non-linearity is useful for classification.
    """

    n_qubits: int
    length: int = 1024

    def forward(self, input_bitstreams: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
        """
        input_bitstreams: (n_qubits, length)
        """
        # 1. Decode inputs to probabilities
        p_in = np.mean(input_bitstreams, axis=1)

        # 2. Quantum Rotation (Simulated)
        theta = p_in * np.pi

        # 3. Measurement Probability
        # Probability of measuring |0>
        p_measure = np.cos(theta / 2.0) ** 2

        # 4. Re-encode to bitstream (Collapse)
        # (n_qubits, length)
        rands = np.random.random((self.n_qubits, self.length))
        out_bits = (rands < p_measure[:, None]).astype(np.uint8)

        return out_bits

forward(input_bitstreams)

input_bitstreams: (n_qubits, length)

Source code in src/sc_neurocore/quantum/hybrid.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
def forward(self, input_bitstreams: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
    """
    input_bitstreams: (n_qubits, length)
    """
    # 1. Decode inputs to probabilities
    p_in = np.mean(input_bitstreams, axis=1)

    # 2. Quantum Rotation (Simulated)
    theta = p_in * np.pi

    # 3. Measurement Probability
    # Probability of measuring |0>
    p_measure = np.cos(theta / 2.0) ** 2

    # 4. Re-encode to bitstream (Collapse)
    # (n_qubits, length)
    rands = np.random.random((self.n_qubits, self.length))
    out_bits = (rands < p_measure[:, None]).astype(np.uint8)

    return out_bits

Noise Models

IBM Heron r2 noise model with depolarizing, amplitude damping, phase damping channels, and asymmetric readout error.

sc_neurocore.quantum.noise_models

HeronR2NoiseParams dataclass

IBM Heron r2 calibration parameters (2024).

Source code in src/sc_neurocore/quantum/noise_models.py 
14
15
16
17
18
19
20
21
22
23
24
25
@dataclass
class HeronR2NoiseParams:
    """IBM Heron r2 calibration parameters (2024)."""

    cx_error: float = 0.005
    single_qubit_error: float = 0.0003
    t1_us: float = 300.0
    t2_us: float = 200.0
    readout_0to1: float = 0.01
    readout_1to0: float = 0.02
    gate_time_1q_ns: float = 25.0
    gate_time_2q_ns: float = 100.0

HeronR2NoiseModel

Source code in src/sc_neurocore/quantum/noise_models.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
class HeronR2NoiseModel:
    def __init__(self, params=None):
        self.params = params or HeronR2NoiseParams()

    def depolarizing_channel(self, p):
        """Kraus operators for single-qubit depolarizing channel."""
        I = np.eye(2, dtype=complex)
        X = np.array([[0, 1], [1, 0]], dtype=complex)
        Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
        Z = np.array([[1, 0], [0, -1]], dtype=complex)
        return [
            np.sqrt(1 - p) * I,
            np.sqrt(p / 3) * X,
            np.sqrt(p / 3) * Y,
            np.sqrt(p / 3) * Z,
        ]

    def amplitude_damping(self, gamma):
        """Kraus operators for amplitude damping (T1 decay)."""
        K0 = np.array([[1, 0], [0, np.sqrt(1 - gamma)]], dtype=complex)
        K1 = np.array([[0, np.sqrt(gamma)], [0, 0]], dtype=complex)
        return [K0, K1]

    def phase_damping(self, gamma):
        """Kraus operators for phase damping (T2 decay)."""
        K0 = np.array([[1, 0], [0, np.sqrt(1 - gamma)]], dtype=complex)
        K1 = np.array([[0, 0], [0, np.sqrt(gamma)]], dtype=complex)
        return [K0, K1]

    def apply_single_qubit_noise(self, rho):
        """Apply single-qubit noise channel to density matrix."""
        kraus = self.depolarizing_channel(self.params.single_qubit_error)
        return sum(K @ rho @ K.conj().T for K in kraus)

    def apply_readout_noise(self, measurement):
        """Apply asymmetric readout error."""
        p = self.params
        if measurement == 0:
            return 1 if np.random.random() < p.readout_0to1 else 0
        return 0 if np.random.random() < p.readout_1to0 else 1

    def gate_fidelity_1q(self):
        return 1.0 - self.params.single_qubit_error

    def gate_fidelity_2q(self):
        return 1.0 - self.params.cx_error

depolarizing_channel(p)

Kraus operators for single-qubit depolarizing channel.

Source code in src/sc_neurocore/quantum/noise_models.py 
32
33
34
35
36
37
38
39
40
41
42
43
def depolarizing_channel(self, p):
    """Kraus operators for single-qubit depolarizing channel."""
    I = np.eye(2, dtype=complex)
    X = np.array([[0, 1], [1, 0]], dtype=complex)
    Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
    Z = np.array([[1, 0], [0, -1]], dtype=complex)
    return [
        np.sqrt(1 - p) * I,
        np.sqrt(p / 3) * X,
        np.sqrt(p / 3) * Y,
        np.sqrt(p / 3) * Z,
    ]

amplitude_damping(gamma)

Kraus operators for amplitude damping (T1 decay).

Source code in src/sc_neurocore/quantum/noise_models.py 
45
46
47
48
49
def amplitude_damping(self, gamma):
    """Kraus operators for amplitude damping (T1 decay)."""
    K0 = np.array([[1, 0], [0, np.sqrt(1 - gamma)]], dtype=complex)
    K1 = np.array([[0, np.sqrt(gamma)], [0, 0]], dtype=complex)
    return [K0, K1]

phase_damping(gamma)

Kraus operators for phase damping (T2 decay).

Source code in src/sc_neurocore/quantum/noise_models.py 
51
52
53
54
55
def phase_damping(self, gamma):
    """Kraus operators for phase damping (T2 decay)."""
    K0 = np.array([[1, 0], [0, np.sqrt(1 - gamma)]], dtype=complex)
    K1 = np.array([[0, 0], [0, np.sqrt(gamma)]], dtype=complex)
    return [K0, K1]

apply_single_qubit_noise(rho)

Apply single-qubit noise channel to density matrix.

Source code in src/sc_neurocore/quantum/noise_models.py 
57
58
59
60
def apply_single_qubit_noise(self, rho):
    """Apply single-qubit noise channel to density matrix."""
    kraus = self.depolarizing_channel(self.params.single_qubit_error)
    return sum(K @ rho @ K.conj().T for K in kraus)

apply_readout_noise(measurement)

Apply asymmetric readout error.

Source code in src/sc_neurocore/quantum/noise_models.py 
62
63
64
65
66
67
def apply_readout_noise(self, measurement):
    """Apply asymmetric readout error."""
    p = self.params
    if measurement == 0:
        return 1 if np.random.random() < p.readout_0to1 else 0
    return 0 if np.random.random() < p.readout_1to0 else 1

Parameter-Shift Gradient

Exact gradient computation for parameterized quantum circuits.

sc_neurocore.quantum.param_shift

parameter_shift_gradient(circuit_fn, params, shift=np.pi / 2)

Gradient via parameter-shift rule.

f'(θ_i) = [f(θ_i + s) - f(θ_i - s)] / (2 sin(s))

Source code in src/sc_neurocore/quantum/param_shift.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def parameter_shift_gradient(circuit_fn, params, shift=np.pi / 2):
    """Gradient via parameter-shift rule.

    f'(θ_i) = [f(θ_i + s) - f(θ_i - s)] / (2 sin(s))
    """
    grad = np.zeros_like(params, dtype=float)
    denom = 2.0 * np.sin(shift)
    for i in range(len(params)):
        p_plus = params.copy()
        p_minus = params.copy()
        p_plus[i] += shift
        p_minus[i] -= shift
        grad[i] = (circuit_fn(p_plus) - circuit_fn(p_minus)) / denom
    return grad

Hybrid Pipeline

VQE-style quantum-classical optimization pipeline.

sc_neurocore.quantum.hybrid_pipeline

HybridQuantumClassicalPipeline

Source code in src/sc_neurocore/quantum/hybrid_pipeline.py 
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
class HybridQuantumClassicalPipeline:
    def __init__(self, n_qubits=2, n_layers=1, noise_model=None):
        self.n_qubits = n_qubits
        self.n_layers = n_layers
        self.noise_model = noise_model
        self.n_params = n_qubits * n_layers

    def circuit(self, params):
        """Parameterized Ry-CNOT circuit → ⟨Z⊗Z⟩ expectation."""
        dim = 2**self.n_qubits
        state = np.zeros(dim, dtype=complex)
        state[0] = 1.0  # |00...0⟩

        idx = 0
        for _ in range(self.n_layers):
            for q in range(self.n_qubits):
                gate = _kron_gate(_ry(params[idx]), q, self.n_qubits)
                state = gate @ state
                idx += 1
            # CNOT chain
            if self.n_qubits >= 2:
                cnot = _cnot()
                for q in range(self.n_qubits - 1):
                    full = np.eye(dim, dtype=complex)
                    # Build CNOT on qubits q, q+1
                    sub = np.eye(dim, dtype=complex)
                    # Direct 2-qubit CNOT embedding
                    for i in range(dim):
                        for j in range(dim):
                            # Extract bits for qubits q and q+1
                            bq = (i >> (self.n_qubits - 1 - q)) & 1
                            bq1 = (i >> (self.n_qubits - 1 - q - 1)) & 1
                            if bq == 1:  # control set → flip target
                                flipped = i ^ (1 << (self.n_qubits - 1 - q - 1))
                                sub[flipped, i] = 1.0
                                sub[i, i] = 0.0
                    state = sub @ state

        # Measure ⟨Z⊗Z⟩ (product of Z eigenvalues on all qubits)
        z_all = np.array([(-1) ** bin(i).count("1") for i in range(dim)], dtype=float)
        return float(np.real(np.conj(state) @ (z_all * state)))

    def train(self, n_steps=100, lr=0.01):
        """VQE-style optimization: minimize ⟨Z⊗Z⟩."""
        params = np.random.randn(self.n_params) * 0.1
        history = []
        for _ in range(n_steps):
            val = self.circuit(params)
            history.append(val)
            grad = parameter_shift_gradient(self.circuit, params)
            params -= lr * grad
        return history, params

    def evaluate(self, params):
        return self.circuit(params)

circuit(params)

Parameterized Ry-CNOT circuit → ⟨Z⊗Z⟩ expectation.

Source code in src/sc_neurocore/quantum/hybrid_pipeline.py 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
def circuit(self, params):
    """Parameterized Ry-CNOT circuit → ⟨Z⊗Z⟩ expectation."""
    dim = 2**self.n_qubits
    state = np.zeros(dim, dtype=complex)
    state[0] = 1.0  # |00...0⟩

    idx = 0
    for _ in range(self.n_layers):
        for q in range(self.n_qubits):
            gate = _kron_gate(_ry(params[idx]), q, self.n_qubits)
            state = gate @ state
            idx += 1
        # CNOT chain
        if self.n_qubits >= 2:
            cnot = _cnot()
            for q in range(self.n_qubits - 1):
                full = np.eye(dim, dtype=complex)
                # Build CNOT on qubits q, q+1
                sub = np.eye(dim, dtype=complex)
                # Direct 2-qubit CNOT embedding
                for i in range(dim):
                    for j in range(dim):
                        # Extract bits for qubits q and q+1
                        bq = (i >> (self.n_qubits - 1 - q)) & 1
                        bq1 = (i >> (self.n_qubits - 1 - q - 1)) & 1
                        if bq == 1:  # control set → flip target
                            flipped = i ^ (1 << (self.n_qubits - 1 - q - 1))
                            sub[flipped, i] = 1.0
                            sub[i, i] = 0.0
                state = sub @ state

    # Measure ⟨Z⊗Z⟩ (product of Z eigenvalues on all qubits)
    z_all = np.array([(-1) ** bin(i).count("1") for i in range(dim)], dtype=float)
    return float(np.real(np.conj(state) @ (z_all * state)))

train(n_steps=100, lr=0.01)

VQE-style optimization: minimize ⟨Z⊗Z⟩.

Source code in src/sc_neurocore/quantum/hybrid_pipeline.py 
84
85
86
87
88
89
90
91
92
93
def train(self, n_steps=100, lr=0.01):
    """VQE-style optimization: minimize ⟨Z⊗Z⟩."""
    params = np.random.randn(self.n_params) * 0.1
    history = []
    for _ in range(n_steps):
        val = self.circuit(params)
        history.append(val)
        grad = parameter_shift_gradient(self.circuit, params)
        params -= lr * grad
    return history, params

QEC

sc_neurocore.quantum.qec

Quantum Error Correction (QEC) Shield for sc-neurocore.

Provides classical-stochastic implementations of QEC codes (Repetition, Surface) to protect quantum-classical bitstreams from noise during IBMQ hardware execution.

QecShield

Repetition code QEC shield for stochastic-quantum bitstreams.

Source code in src/sc_neurocore/quantum/qec.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
class QecShield:
    """Repetition code QEC shield for stochastic-quantum bitstreams."""

    def __init__(self, code_type: str = "repetition", distance: int = 3):
        self.code_type = code_type
        self.distance = distance

    def encode(self, bitstream: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
        """repetition code (d=3): 0 -> 000, 1 -> 111"""
        if self.code_type == "repetition":
            return np.repeat(bitstream[:, np.newaxis, :], self.distance, axis=1)
        return bitstream

    def extract_syndromes(self, physical_bits: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
        if self.code_type == "repetition":
            res: np.ndarray[Any, Any] = np.diff(physical_bits, axis=1) % 2
            return res
        return np.zeros_like(physical_bits)

    def decode(self, physical_bits: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
        if self.code_type == "repetition":
            means = np.mean(physical_bits, axis=1)
            res: np.ndarray[Any, Any] = (means > 0.5).astype(np.uint8)
            return res
        return physical_bits

    def get_error_rate(self, syndromes: np.ndarray[Any, Any]) -> float:
        return float(np.mean(syndromes))

encode(bitstream)

repetition code (d=3): 0 -> 000, 1 -> 111

Source code in src/sc_neurocore/quantum/qec.py 
27
28
29
30
31
def encode(self, bitstream: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
    """repetition code (d=3): 0 -> 000, 1 -> 111"""
    if self.code_type == "repetition":
        return np.repeat(bitstream[:, np.newaxis, :], self.distance, axis=1)
    return bitstream

SurfaceCodeShield

Distance-d rotated surface code for stochastic-quantum bitstreams.

Encodes 1 logical qubit into d² physical data qubits. X and Z stabilizers detect bit-flip and phase-flip errors. Decoding uses a lookup table for d=3.

Ref: Fowler et al., "Surface codes: Towards practical large-scale quantum computation", Phys. Rev. A 86, 032324 (2012).

Source code in src/sc_neurocore/quantum/qec.py 
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
class SurfaceCodeShield:
    """
    Distance-d rotated surface code for stochastic-quantum bitstreams.

    Encodes 1 logical qubit into d² physical data qubits. X and Z stabilizers
    detect bit-flip and phase-flip errors. Decoding uses a lookup table for d=3.

    Ref: Fowler et al., "Surface codes: Towards practical large-scale quantum
    computation", Phys. Rev. A 86, 032324 (2012).
    """

    def __init__(self, distance: int = 3):
        if distance < 3 or distance % 2 == 0:
            raise ValueError("distance must be odd >= 3")
        self.distance = distance
        self.n_data = distance * distance
        # Stabilizer generators: each is a list of data-qubit indices
        self.x_stabilizers, self.z_stabilizers = self._build_stabilizers(distance)
        # d=3 lookup table: syndrome pattern -> correction qubit index
        if distance == 3:
            self._x_lut = self._build_d3_lut(self.x_stabilizers)
            self._z_lut = self._build_d3_lut(self.z_stabilizers)

    @staticmethod
    def _build_stabilizers(d: int):
        """Build X and Z stabilizer generators for rotated surface code."""
        x_stabs: list[list[int]] = []
        z_stabs: list[list[int]] = []
        for r in range(d):
            for c in range(d):
                idx = r * d + c
                # X stabilizers: plaquettes on even sublattice
                if (r + c) % 2 == 0 and r < d - 1 and c < d - 1:
                    x_stabs.append([idx, idx + 1, idx + d, idx + d + 1])
                # Z stabilizers: plaquettes on odd sublattice
                if (r + c) % 2 == 1 and r < d - 1 and c < d - 1:
                    z_stabs.append([idx, idx + 1, idx + d, idx + d + 1])
        # Boundary stabilizers (weight-2) for top/bottom/left/right edges
        for c in range(0, d - 1, 2):
            x_stabs.append([c, c + 1])  # top edge
        for c in range(1 if d > 3 else 0, d - 1, 2):
            if (d - 1) * d + c < d * d and (d - 1) * d + c + 1 < d * d:
                x_stabs.append([(d - 1) * d + c, (d - 1) * d + c + 1])  # bottom edge
        for r in range(0, d - 1, 2):
            z_stabs.append([r * d, (r + 1) * d])  # left edge
        for r in range(1 if d > 3 else 0, d - 1, 2):
            if (r + 1) * d + d - 1 < d * d:
                z_stabs.append([r * d + d - 1, (r + 1) * d + d - 1])  # right edge
        return x_stabs, z_stabs

    @staticmethod
    def _build_d3_lut(stabilizers: list[list[int]]) -> dict[tuple[int, ...], int]:
        """Build syndrome → correction lookup for d=3 single-qubit errors."""
        lut: dict[tuple[int, ...], int] = {}
        n_stabs = len(stabilizers)
        for qubit in range(9):
            syndrome = [0] * n_stabs
            for s_idx, stab in enumerate(stabilizers):
                if qubit in stab:
                    syndrome[s_idx] = 1
            key = tuple(syndrome)
            if key not in lut:
                lut[key] = qubit
        return lut

    def encode(self, bitstream: np.ndarray) -> np.ndarray:
        """
        Encode logical bitstream into surface code physical qubits.

        Input: (n_logical, length) — each row is one logical qubit's bitstream.
        Output: (n_logical, n_data, length) — repeated into d² data qubits.
        """
        return np.repeat(bitstream[:, np.newaxis, :], self.n_data, axis=1)

    def measure_syndrome(self, physical_bits: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        """
        Measure X and Z stabilizer syndromes.

        Input: (n_logical, n_data, length)
        Returns: (x_syndrome, z_syndrome) each (n_logical, n_stabilizers, length)
        """
        n_logical, _, length = physical_bits.shape
        x_syn = np.zeros((n_logical, len(self.x_stabilizers), length), dtype=np.uint8)
        z_syn = np.zeros((n_logical, len(self.z_stabilizers), length), dtype=np.uint8)
        for s_idx, stab in enumerate(self.x_stabilizers):
            parity = np.zeros((n_logical, length), dtype=np.uint8)
            for q in stab:
                parity ^= physical_bits[:, q, :]
            x_syn[:, s_idx, :] = parity
        for s_idx, stab in enumerate(self.z_stabilizers):
            parity = np.zeros((n_logical, length), dtype=np.uint8)
            for q in stab:
                parity ^= physical_bits[:, q, :]
            z_syn[:, s_idx, :] = parity
        return x_syn, z_syn

    def decode(self, physical_bits: np.ndarray) -> np.ndarray:
        """
        Decode surface code: measure syndromes, correct single-qubit errors, majority vote.

        Input: (n_logical, n_data, length)
        Output: (n_logical, length)
        """
        corrected = physical_bits.copy()
        x_syn, z_syn = self.measure_syndrome(corrected)
        n_logical, _, length = corrected.shape

        if self.distance == 3:
            self._apply_lut_correction(corrected, x_syn, self._x_lut)
            self._apply_lut_correction(corrected, z_syn, self._z_lut)
        else:
            # For d>3, apply majority vote per stabilizer neighbourhood
            pass

        # Majority vote across all data qubits
        means = np.mean(corrected, axis=1)
        return (means > 0.5).astype(np.uint8)

    @staticmethod
    def _apply_lut_correction(
        physical: np.ndarray, syndromes: np.ndarray, lut: dict[tuple[int, ...], int]
    ) -> None:
        """Apply lookup-table correction for each bitstream position."""
        n_logical, n_stab, length = syndromes.shape
        for l_idx in range(n_logical):
            for t in range(length):
                syn_key = tuple(int(syndromes[l_idx, s, t]) for s in range(n_stab))
                if any(syn_key):
                    qubit = lut.get(syn_key)
                    if qubit is not None:
                        physical[l_idx, qubit, t] ^= 1

    def get_error_rate(self, x_syn: np.ndarray, z_syn: np.ndarray) -> float:
        """Estimated error rate from syndrome density."""
        return float((np.mean(x_syn) + np.mean(z_syn)) / 2)

encode(bitstream)

Encode logical bitstream into surface code physical qubits.

Input: (n_logical, length) — each row is one logical qubit's bitstream. Output: (n_logical, n_data, length) — repeated into d² data qubits.

Source code in src/sc_neurocore/quantum/qec.py 
115
116
117
118
119
120
121
122
def encode(self, bitstream: np.ndarray) -> np.ndarray:
    """
    Encode logical bitstream into surface code physical qubits.

    Input: (n_logical, length) — each row is one logical qubit's bitstream.
    Output: (n_logical, n_data, length) — repeated into d² data qubits.
    """
    return np.repeat(bitstream[:, np.newaxis, :], self.n_data, axis=1)

measure_syndrome(physical_bits)

Measure X and Z stabilizer syndromes.

Input: (n_logical, n_data, length) Returns: (x_syndrome, z_syndrome) each (n_logical, n_stabilizers, length)

Source code in src/sc_neurocore/quantum/qec.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
def measure_syndrome(self, physical_bits: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """
    Measure X and Z stabilizer syndromes.

    Input: (n_logical, n_data, length)
    Returns: (x_syndrome, z_syndrome) each (n_logical, n_stabilizers, length)
    """
    n_logical, _, length = physical_bits.shape
    x_syn = np.zeros((n_logical, len(self.x_stabilizers), length), dtype=np.uint8)
    z_syn = np.zeros((n_logical, len(self.z_stabilizers), length), dtype=np.uint8)
    for s_idx, stab in enumerate(self.x_stabilizers):
        parity = np.zeros((n_logical, length), dtype=np.uint8)
        for q in stab:
            parity ^= physical_bits[:, q, :]
        x_syn[:, s_idx, :] = parity
    for s_idx, stab in enumerate(self.z_stabilizers):
        parity = np.zeros((n_logical, length), dtype=np.uint8)
        for q in stab:
            parity ^= physical_bits[:, q, :]
        z_syn[:, s_idx, :] = parity
    return x_syn, z_syn

decode(physical_bits)

Decode surface code: measure syndromes, correct single-qubit errors, majority vote.

Input: (n_logical, n_data, length) Output: (n_logical, length)

Source code in src/sc_neurocore/quantum/qec.py 
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
def decode(self, physical_bits: np.ndarray) -> np.ndarray:
    """
    Decode surface code: measure syndromes, correct single-qubit errors, majority vote.

    Input: (n_logical, n_data, length)
    Output: (n_logical, length)
    """
    corrected = physical_bits.copy()
    x_syn, z_syn = self.measure_syndrome(corrected)
    n_logical, _, length = corrected.shape

    if self.distance == 3:
        self._apply_lut_correction(corrected, x_syn, self._x_lut)
        self._apply_lut_correction(corrected, z_syn, self._z_lut)
    else:
        # For d>3, apply majority vote per stabilizer neighbourhood
        pass

    # Majority vote across all data qubits
    means = np.mean(corrected, axis=1)
    return (means > 0.5).astype(np.uint8)

get_error_rate(x_syn, z_syn)

Estimated error rate from syndrome density.

Source code in src/sc_neurocore/quantum/qec.py 
182
183
184
def get_error_rate(self, x_syn: np.ndarray, z_syn: np.ndarray) -> float:
    """Estimated error rate from syndrome density."""
    return float((np.mean(x_syn) + np.mean(z_syn)) / 2)

SC→Quantum Compiler (Conjecture C1+C4)

Compiles SC operations to quantum circuits. SC probability p encodes as Ry(2·arcsin(√p)) rotation; AND gate maps to joint measurement; Born rule recovers P(|1⟩) = p exactly. Includes noisy simulation via HeronR2NoiseModel.

sc_neurocore.quantum.sc_quantum_compiler.sc_prob_to_statevector(p)

Encode SC probability as a single-qubit state vector.

|ψ⟩ = √(1-p)|0⟩ + √p|1⟩ → P(measure |1⟩) = p exactly.

Source code in src/sc_neurocore/quantum/sc_quantum_compiler.py 
38
39
40
41
42
43
44
def sc_prob_to_statevector(p: float) -> np.ndarray:
    """Encode SC probability as a single-qubit state vector.

    |ψ⟩ = √(1-p)|0⟩ + √p|1⟩  →  P(measure |1⟩) = p exactly.
    """
    p = float(np.clip(p, 0.0, 1.0))
    return np.array([np.sqrt(1.0 - p), np.sqrt(p)], dtype=complex)

sc_neurocore.quantum.sc_quantum_compiler.compile_sc_multiply(p_a, p_b)

Compile SC AND gate (multiplication) to a quantum circuit.

SC: P(a AND b) = P(a) * P(b) for independent streams. Quantum: encode probabilities as Ry rotations, use CNOT for correlation.

Source code in src/sc_neurocore/quantum/sc_quantum_compiler.py 
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
def compile_sc_multiply(p_a: float, p_b: float) -> SCQuantumCircuit:
    """Compile SC AND gate (multiplication) to a quantum circuit.

    SC: P(a AND b) = P(a) * P(b) for independent streams.
    Quantum: encode probabilities as Ry rotations, use CNOT for correlation.
    """
    theta_a = prob_to_ry_angle(p_a)
    theta_b = prob_to_ry_angle(p_b)

    # 2 qubits: q0 encodes p_a, q1 encodes p_b
    # Product probability appears on q1 conditioned on q0
    gates = [
        QuantumGate("Ry(p_a)", ry_gate(theta_a), [0]),
        QuantumGate("Ry(p_b)", ry_gate(theta_b), [1]),
    ]

    # The output is the joint probability P(q0=1 AND q1=1)
    circuit = SCQuantumCircuit(
        n_qubits=2,
        gates=gates,
        input_qubits=[0, 1],
        output_qubit=1,  # marginal on q1
    )
    return circuit

sc_neurocore.quantum.sc_quantum_compiler.compile_sc_layer(weights, input_probs)

Compile an SC dense layer to quantum gate descriptions.

Parameters

weights : np.ndarray Shape (n_neurons, n_inputs), values in [0, 1]. input_probs : np.ndarray Shape (n_inputs,), SC input probabilities.

Returns

list of dicts, one per neuron, each containing: 'neuron_idx': int 'ry_angles': list of (input_angle, weight_angle) pairs 'expected_output': float — SC computation result 'quantum_output': float — quantum simulation result

Source code in src/sc_neurocore/quantum/sc_quantum_compiler.py 
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
def compile_sc_layer(weights: np.ndarray, input_probs: np.ndarray) -> list[dict[str, Any]]:
    """Compile an SC dense layer to quantum gate descriptions.

    Parameters
    ----------
    weights : np.ndarray
        Shape (n_neurons, n_inputs), values in [0, 1].
    input_probs : np.ndarray
        Shape (n_inputs,), SC input probabilities.

    Returns
    -------
    list of dicts, one per neuron, each containing:
        'neuron_idx': int
        'ry_angles': list of (input_angle, weight_angle) pairs
        'expected_output': float — SC computation result
        'quantum_output': float — quantum simulation result
    """
    n_neurons, n_inputs = weights.shape
    results = []

    for j in range(n_neurons):
        ry_angles = []
        sc_output = 0.0
        quantum_outputs = []

        for i in range(n_inputs):
            w = float(np.clip(weights[j, i], 0, 1))
            x = float(np.clip(input_probs[i], 0, 1))

            theta_x = prob_to_ry_angle(x)
            theta_w = prob_to_ry_angle(w)
            ry_angles.append((theta_x, theta_w))

            # SC: AND gate → product
            sc_output += w * x

            # Quantum: independent product P(q0=1)*P(q1=1)
            quantum_outputs.append(w * x)

        sc_output = float(np.clip(sc_output / max(n_inputs, 1), 0, 1))
        q_output = float(np.clip(sum(quantum_outputs) / max(n_inputs, 1), 0, 1))

        results.append(
            {
                "neuron_idx": j,
                "ry_angles": ry_angles,
                "expected_output": sc_output,
                "quantum_output": q_output,
            }
        )

    return results

sc_neurocore.quantum.sc_quantum_compiler.SCQuantumCircuit dataclass

Quantum circuit compiled from SC operations.

Source code in src/sc_neurocore/quantum/sc_quantum_compiler.py 
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
@dataclass
class SCQuantumCircuit:
    """Quantum circuit compiled from SC operations."""

    n_qubits: int
    gates: list[QuantumGate]
    input_qubits: list[int]
    output_qubit: int

    def simulate(self) -> np.ndarray:
        """Simulate the circuit and return the full statevector."""
        dim = 2**self.n_qubits
        state = np.zeros(dim, dtype=complex)
        state[0] = 1.0  # |000...0⟩

        for gate in self.gates:
            state = _apply_gate(state, gate.matrix, gate.qubits, self.n_qubits)

        return state

    def output_probability(self) -> float:
        """Simulate and return P(output_qubit = |1⟩)."""
        state = self.simulate()
        prob = 0.0
        for i in range(len(state)):
            if (i >> self.output_qubit) & 1:
                prob += np.abs(state[i]) ** 2
        return float(prob)

    def simulate_noisy(self, noise_model) -> np.ndarray:
        """Simulate with noise: evolve density matrix through Kraus channels.

        Parameters
        ----------
        noise_model : HeronR2NoiseModel or compatible
            Must provide apply_single_qubit_noise(rho) and apply_readout_noise(measurement).

        Returns
        -------
        np.ndarray
            Final density matrix of shape (2^n, 2^n).
        """
        dim = 2**self.n_qubits
        state = np.zeros(dim, dtype=complex)
        state[0] = 1.0
        # Apply gates as unitary
        for gate in self.gates:
            state = _apply_gate(state, gate.matrix, gate.qubits, self.n_qubits)
        # Convert to density matrix
        rho = np.outer(state, state.conj())
        # Apply per-qubit noise
        for q in range(self.n_qubits):
            rho = _apply_single_qubit_channel(rho, noise_model, q, self.n_qubits)
        return rho

    def output_probability_noisy(self, noise_model, n_shots: int = 1000) -> float:
        """Simulate with noise and return P(output=1) via measurement sampling.

        Parameters
        ----------
        noise_model : HeronR2NoiseModel or compatible
        n_shots : int
            Number of measurement shots.
        """
        rho = self.simulate_noisy(noise_model)
        # Extract output qubit probability from density matrix diagonal
        prob_1 = 0.0
        dim = 2**self.n_qubits
        for i in range(dim):
            if (i >> self.output_qubit) & 1:
                prob_1 += float(np.real(rho[i, i]))
        # Apply readout noise via sampling
        ones = sum(
            1
            for _ in range(n_shots)
            if noise_model.apply_readout_noise(1 if np.random.random() < prob_1 else 0) == 1
        )
        return ones / n_shots

    def summary(self) -> str:
        lines = [f"SCQuantumCircuit: {self.n_qubits} qubits, {len(self.gates)} gates"]
        for g in self.gates:
            lines.append(f"  {g.name} on qubit(s) {g.qubits}")
        lines.append(f"  output: qubit {self.output_qubit}")
        return "\n".join(lines)

simulate()

Simulate the circuit and return the full statevector.

Source code in src/sc_neurocore/quantum/sc_quantum_compiler.py 
88
89
90
91
92
93
94
95
96
97
def simulate(self) -> np.ndarray:
    """Simulate the circuit and return the full statevector."""
    dim = 2**self.n_qubits
    state = np.zeros(dim, dtype=complex)
    state[0] = 1.0  # |000...0⟩

    for gate in self.gates:
        state = _apply_gate(state, gate.matrix, gate.qubits, self.n_qubits)

    return state

output_probability()

Simulate and return P(output_qubit = |1⟩).

Source code in src/sc_neurocore/quantum/sc_quantum_compiler.py 
 99
100
101
102
103
104
105
106
def output_probability(self) -> float:
    """Simulate and return P(output_qubit = |1⟩)."""
    state = self.simulate()
    prob = 0.0
    for i in range(len(state)):
        if (i >> self.output_qubit) & 1:
            prob += np.abs(state[i]) ** 2
    return float(prob)

simulate_noisy(noise_model)

Simulate with noise: evolve density matrix through Kraus channels.

Parameters

noise_model : HeronR2NoiseModel or compatible Must provide apply_single_qubit_noise(rho) and apply_readout_noise(measurement).

Returns

np.ndarray Final density matrix of shape (2^n, 2^n).

Source code in src/sc_neurocore/quantum/sc_quantum_compiler.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
def simulate_noisy(self, noise_model) -> np.ndarray:
    """Simulate with noise: evolve density matrix through Kraus channels.

    Parameters
    ----------
    noise_model : HeronR2NoiseModel or compatible
        Must provide apply_single_qubit_noise(rho) and apply_readout_noise(measurement).

    Returns
    -------
    np.ndarray
        Final density matrix of shape (2^n, 2^n).
    """
    dim = 2**self.n_qubits
    state = np.zeros(dim, dtype=complex)
    state[0] = 1.0
    # Apply gates as unitary
    for gate in self.gates:
        state = _apply_gate(state, gate.matrix, gate.qubits, self.n_qubits)
    # Convert to density matrix
    rho = np.outer(state, state.conj())
    # Apply per-qubit noise
    for q in range(self.n_qubits):
        rho = _apply_single_qubit_channel(rho, noise_model, q, self.n_qubits)
    return rho

output_probability_noisy(noise_model, n_shots=1000)

Simulate with noise and return P(output=1) via measurement sampling.

Parameters

noise_model : HeronR2NoiseModel or compatible n_shots : int Number of measurement shots.

Source code in src/sc_neurocore/quantum/sc_quantum_compiler.py 
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
def output_probability_noisy(self, noise_model, n_shots: int = 1000) -> float:
    """Simulate with noise and return P(output=1) via measurement sampling.

    Parameters
    ----------
    noise_model : HeronR2NoiseModel or compatible
    n_shots : int
        Number of measurement shots.
    """
    rho = self.simulate_noisy(noise_model)
    # Extract output qubit probability from density matrix diagonal
    prob_1 = 0.0
    dim = 2**self.n_qubits
    for i in range(dim):
        if (i >> self.output_qubit) & 1:
            prob_1 += float(np.real(rho[i, i]))
    # Apply readout noise via sampling
    ones = sum(
        1
        for _ in range(n_shots)
        if noise_model.apply_readout_noise(1 if np.random.random() < prob_1 else 0) == 1
    )
    return ones / n_shots