From ec0142203770436cc7e4aca5f01ca2d67cebcc46 Mon Sep 17 00:00:00 2001 From: BrunoLiegiBastonLiegi Date: Tue, 25 Nov 2025 11:36:56 +0100 Subject: [PATCH 1/4] feat: improved measurements sampling --- src/qibo/backends/_clifford_operations.py | 29 ++++++++++++++++++++--- src/qibo/backends/clifford.py | 14 +++-------- tests/test_backends_clifford.py | 2 +- 3 files changed, 30 insertions(+), 15 deletions(-) diff --git a/src/qibo/backends/_clifford_operations.py b/src/qibo/backends/_clifford_operations.py index e893631188..3f79db813d 100644 --- a/src/qibo/backends/_clifford_operations.py +++ b/src/qibo/backends/_clifford_operations.py @@ -3,6 +3,8 @@ import numpy as np from scipy import sparse +from qibo.config import raise_error + name = "numpy" @@ -494,8 +496,7 @@ def _init_state_for_measurements(state, nqubits, collapse): return state.copy() -# valid for a standard basis measurement only -def M(state, qubits, nqubits, collapse=False): +def _M(state, qubits, nqubits, collapse=False): sample = [] state = _init_state_for_measurements(state, nqubits, collapse) # TODO: parallelize this and get rid of the loop @@ -503,16 +504,38 @@ def M(state, qubits, nqubits, collapse=False): p = state[nqubits:-1, q].nonzero()[0] # random outcome, affects the state if len(p) > 0: - state, outcome = _random_outcome(state, p, q, nqubits) + state, _ = _random_outcome(state, p, q, nqubits) + outcome = None # determined outcome, state unchanged else: _, outcome = _determined_outcome(state, q, nqubits) + outcome = int(outcome) sample.append(outcome) if collapse: state = _packbits(state, axis=0) return sample +# valid for a standard basis measurement only +def M(state, qubits, nqubits, collapse=False, nshots=1): + if collapse and nshots != 1: + raise_error( + RuntimeError, + "Cannot generate multiple shots with a collapsing measurement.", + ) + sample = _M(state, qubits, nqubits, collapse) + samples = [ + ( + np.random.randint(2, size=nshots).tolist() + if outcome is None + else nshots * [outcome] + ) + for outcome in sample + ] + samples = list(zip(*samples)) + return samples + + def cast(x, dtype=None, copy: bool = False): if dtype is None: dtype = "complex128" diff --git a/src/qibo/backends/clifford.py b/src/qibo/backends/clifford.py index 237267ed41..b9a31507e2 100644 --- a/src/qibo/backends/clifford.py +++ b/src/qibo/backends/clifford.py @@ -410,14 +410,6 @@ def execute_circuit( # pylint: disable=R1710 if self.platform == "stim": return self._execute_circuit_stim(circuit, initial_state, nshots) - for gate in circuit.queue: - if ( - not gate.clifford - and not gate.__class__.__name__ == "M" - and not isinstance(gate, gates.PauliNoiseChannel) - ): - raise_error(RuntimeError, "Circuit contains non-Clifford gates.") - if circuit.repeated_execution and nshots != 1: return self.execute_circuit_repeated(circuit, nshots, initial_state) @@ -531,10 +523,10 @@ def sample_shots( qubits = tuple(qubits) if collapse: - samples = [self.engine.M(state, qubits, nqubits) for _ in range(nshots - 1)] - samples.append(self.engine.M(state, qubits, nqubits, collapse)) + samples = self.engine.M(state, qubits, nqubits, nshots=nshots - 1) + samples.append(self.engine.M(state, qubits, nqubits, collapse, nshots=1)) else: - samples = [self.engine.M(state, qubits, nqubits) for _ in range(nshots)] + samples = self.engine.M(state, qubits, nqubits, nshots=nshots) return self.engine.cast(samples, dtype=int) def symplectic_matrix_to_generators( diff --git a/tests/test_backends_clifford.py b/tests/test_backends_clifford.py index d5cb19fc5b..d333849f78 100644 --- a/tests/test_backends_clifford.py +++ b/tests/test_backends_clifford.py @@ -287,7 +287,7 @@ def test_non_clifford_error(backend): clifford_bkd = construct_clifford_backend(backend) c = Circuit(1) c.add(gates.T(0)) - with pytest.raises(RuntimeError) as excinfo: + with pytest.raises(AttributeError) as excinfo: clifford_bkd.execute_circuit(c) assert str(excinfo.value) == "Circuit contains non-Clifford gates." From acb93dbb01b4b2ce0808089e7e02328c68adb558 Mon Sep 17 00:00:00 2001 From: BrunoLiegiBastonLiegi Date: Tue, 25 Nov 2025 12:16:25 +0100 Subject: [PATCH 2/4] not working actually --- src/qibo/backends/#_clifford_operations.py# | 585 ++++++++++++++++++++ src/qibo/backends/.#_clifford_operations.py | 1 + src/qibo/backends/_clifford_operations.py | 16 +- tests/test_backends_clifford.py | 3 +- 4 files changed, 596 insertions(+), 9 deletions(-) create mode 100644 src/qibo/backends/#_clifford_operations.py# create mode 120000 src/qibo/backends/.#_clifford_operations.py diff --git a/src/qibo/backends/#_clifford_operations.py# b/src/qibo/backends/#_clifford_operations.py# new file mode 100644 index 0000000000..786c2e7e4e --- /dev/null +++ b/src/qibo/backends/#_clifford_operations.py# @@ -0,0 +1,585 @@ +from functools import cache, reduce + +import numpy as np +from scipy import sparse + +from qibo.config import raise_error + +name = "numpy" + + +def _get_rxz(symplectic_matrix, nqubits): + return ( + symplectic_matrix[:, -1], + symplectic_matrix[:, :nqubits], + symplectic_matrix[:, nqubits:-1], + ) + + +def I(symplectic_matrix, q, nqubits): + return symplectic_matrix + + +def H(symplectic_matrix, q, nqubits): + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = r ^ (x[:, q] & z[:, q]) + symplectic_matrix[:, [q, nqubits + q]] = symplectic_matrix[:, [nqubits + q, q]] + return symplectic_matrix + + +def CNOT(symplectic_matrix, control_q, target_q, nqubits): + ind_zt = nqubits + target_q + ind_zc = nqubits + control_q + r = symplectic_matrix[:, -1] + xcq = symplectic_matrix[:, control_q] + xtq = symplectic_matrix[:, target_q] + ztq = symplectic_matrix[:, ind_zt] + zcq = symplectic_matrix[:, ind_zc] + symplectic_matrix[:, -1] = r ^ (xcq & ztq) & (xtq ^ ~zcq) + symplectic_matrix[:, target_q] = xtq ^ xcq + symplectic_matrix[:, ind_zc] = zcq ^ ztq + return symplectic_matrix + + +def CZ(symplectic_matrix, control_q, target_q, nqubits): + """Decomposition --> H-CNOT-H""" + ind_zt = nqubits + target_q + ind_zc = nqubits + control_q + r = symplectic_matrix[:, -1] + xcq = symplectic_matrix[:, control_q] + xtq = symplectic_matrix[:, target_q] + ztq = symplectic_matrix[:, ind_zt] + zcq = symplectic_matrix[:, ind_zc] + ztq_xor_xcq = ztq ^ xcq + symplectic_matrix[:, -1] = ( + r ^ (xtq & ztq) ^ (xcq & xtq & (ztq ^ ~zcq)) ^ (xtq & ztq_xor_xcq) + ) + z_control_q = xtq ^ zcq + z_target_q = ztq_xor_xcq + symplectic_matrix[:, ind_zc] = z_control_q + symplectic_matrix[:, ind_zt] = z_target_q + return symplectic_matrix + + +def S(symplectic_matrix, q, nqubits): + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = r ^ (x[:, q] & z[:, q]) + symplectic_matrix[:, q + nqubits] = z[:, q] ^ x[:, q] + return symplectic_matrix + + +def Z(symplectic_matrix, q, nqubits): + """Decomposition --> S-S""" + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = r ^ ((x[:, q] & z[:, q]) ^ x[:, q] & (z[:, q] ^ x[:, q])) + return symplectic_matrix + + +def X(symplectic_matrix, q, nqubits): + """Decomposition --> H-S-S-H""" + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = r ^ (z[:, q] & (z[:, q] ^ x[:, q])) ^ (z[:, q] & x[:, q]) + return symplectic_matrix + + +def Y(symplectic_matrix, q, nqubits): + """Decomposition --> S-S-H-S-S-H""" + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = ( + r ^ (z[:, q] & (z[:, q] ^ x[:, q])) ^ (x[:, q] & (z[:, q] ^ x[:, q])) + ) + return symplectic_matrix + + +def SX(symplectic_matrix, q, nqubits): + """Decomposition --> H-S-H""" + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = r ^ (z[:, q] & (z[:, q] ^ x[:, q])) + symplectic_matrix[:, q] = z[:, q] ^ x[:, q] + return symplectic_matrix + + +def SDG(symplectic_matrix, q, nqubits): + """Decomposition --> S-S-S""" + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = r ^ (x[:, q] & (z[:, q] ^ x[:, q])) + symplectic_matrix[:, nqubits + q] = z[:, q] ^ x[:, q] + return symplectic_matrix + + +def SXDG(symplectic_matrix, q, nqubits): + """Decomposition --> H-S-S-S-H""" + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = r ^ (z[:, q] & x[:, q]) + symplectic_matrix[:, q] = z[:, q] ^ x[:, q] + return symplectic_matrix + + +def RX(symplectic_matrix, q, nqubits, theta): + if theta % (2 * np.pi) == 0: + return I(symplectic_matrix, q, nqubits) + elif (theta / np.pi - 1) % 2 == 0: + return X(symplectic_matrix, q, nqubits) + elif (theta / (np.pi / 2) - 1) % 4 == 0: + return SX(symplectic_matrix, q, nqubits) + else: # theta == 3*pi/2 + 2*n*pi + return SXDG(symplectic_matrix, q, nqubits) + + +def RZ(symplectic_matrix, q, nqubits, theta): + if theta % (2 * np.pi) == 0: + return I(symplectic_matrix, q, nqubits) + elif (theta / np.pi - 1) % 2 == 0: + return Z(symplectic_matrix, q, nqubits) + elif (theta / (np.pi / 2) - 1) % 4 == 0: + return S(symplectic_matrix, q, nqubits) + else: # theta == 3*pi/2 + 2*n*pi + return SDG(symplectic_matrix, q, nqubits) + + +def RY_pi(symplectic_matrix, q, nqubits): + """Decomposition --> H-S-S""" + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = r ^ (x[:, q] & (z[:, q] ^ x[:, q])) + symplectic_matrix[:, [nqubits + q, q]] = symplectic_matrix[:, [q, nqubits + q]] + return symplectic_matrix + + +def RY_3pi_2(symplectic_matrix, q, nqubits): + """Decomposition --> H-S-S-H-S-S-H-S-S""" + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = r ^ (z[:, q] & (z[:, q] ^ x[:, q])) + symplectic_matrix[:, [nqubits + q, q]] = symplectic_matrix[:, [q, nqubits + q]] + return symplectic_matrix + + +def RY(symplectic_matrix, q, nqubits, theta): + if theta % (2 * np.pi) == 0: + return I(symplectic_matrix, q, nqubits) + elif (theta / np.pi - 1) % 2 == 0: + return Y(symplectic_matrix, q, nqubits) + elif (theta / (np.pi / 2) - 1) % 4 == 0: + """Decomposition --> H-S-S""" + return RY_pi(symplectic_matrix, q, nqubits) + else: # theta == 3*pi/2 + 2*n*pi + """Decomposition --> H-S-S-H-S-S-H-S-S""" + return RY_3pi_2(symplectic_matrix, q, nqubits) + + +def GPI2(symplectic_matrix, q, nqubits, phi): + if phi % (2 * np.pi) == 0: + return RX(symplectic_matrix, q, nqubits, np.pi / 2) + + if (phi / np.pi - 1) % 2 == 0: + return RX(symplectic_matrix, q, nqubits, -np.pi / 2) + + if (phi / (np.pi / 2) - 1) % 4 == 0: + return RY(symplectic_matrix, q, nqubits, np.pi / 2) + + # theta == 3*pi/2 + 2*n*pi + return RY(symplectic_matrix, q, nqubits, -np.pi / 2) + + +def SWAP(symplectic_matrix, control_q, target_q, nqubits): + """Decomposition --> CNOT-CNOT-CNOT""" + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = ( + r + ^ (x[:, control_q] & z[:, target_q] & (x[:, target_q] ^ ~z[:, control_q])) + ^ ( + (x[:, target_q] ^ x[:, control_q]) + & (z[:, target_q] ^ z[:, control_q]) + & (z[:, target_q] ^ ~x[:, control_q]) + ) + ^ ( + x[:, target_q] + & z[:, control_q] + & (x[:, control_q] ^ x[:, target_q] ^ z[:, control_q] ^ ~z[:, target_q]) + ) + ) + symplectic_matrix[ + :, [control_q, target_q, nqubits + control_q, nqubits + target_q] + ] = symplectic_matrix[ + :, [target_q, control_q, nqubits + target_q, nqubits + control_q] + ] + return symplectic_matrix + + +def iSWAP(symplectic_matrix, control_q, target_q, nqubits): + """Decomposition --> H-CNOT-CNOT-H-S-S""" + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = ( + r + ^ (x[:, target_q] & z[:, target_q]) + ^ (x[:, control_q] & z[:, control_q]) + ^ (x[:, control_q] & (z[:, control_q] ^ x[:, control_q])) + ^ ( + (z[:, control_q] ^ x[:, control_q]) + & (z[:, target_q] ^ x[:, target_q]) + & (x[:, target_q] ^ ~x[:, control_q]) + ) + ^ ( + (x[:, target_q] ^ z[:, control_q] ^ x[:, control_q]) + & (x[:, target_q] ^ z[:, target_q] ^ x[:, control_q]) + & (x[:, target_q] ^ z[:, target_q] ^ x[:, control_q] ^ ~z[:, control_q]) + ) + ^ (x[:, control_q] & (x[:, target_q] ^ x[:, control_q] ^ z[:, control_q])) + ) + z_control_q = x[:, target_q] ^ z[:, target_q] ^ x[:, control_q] + z_target_q = x[:, target_q] ^ z[:, control_q] ^ x[:, control_q] + symplectic_matrix[:, nqubits + control_q] = z_control_q + symplectic_matrix[:, nqubits + target_q] = z_target_q + symplectic_matrix[:, [control_q, target_q]] = symplectic_matrix[ + :, [target_q, control_q] + ] + return symplectic_matrix + + +def FSWAP(symplectic_matrix, control_q, target_q, nqubits): + """Decomposition --> X-CNOT-RY-CNOT-RY-CNOT-CNOT-X""" + symplectic_matrix = X(symplectic_matrix, target_q, nqubits) + symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) + symplectic_matrix = RY(symplectic_matrix, control_q, nqubits, np.pi / 2) + symplectic_matrix = CNOT(symplectic_matrix, target_q, control_q, nqubits) + symplectic_matrix = RY(symplectic_matrix, control_q, nqubits, -np.pi / 2) + symplectic_matrix = CNOT(symplectic_matrix, target_q, control_q, nqubits) + symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) + return X(symplectic_matrix, control_q, nqubits) + + +def CY(symplectic_matrix, control_q, target_q, nqubits): + """Decomposition --> S-CNOT-SDG""" + r, x, z = _get_rxz(symplectic_matrix, nqubits) + symplectic_matrix[:, -1] = ( + r + ^ (x[:, target_q] & (z[:, target_q] ^ x[:, target_q])) + ^ ( + x[:, control_q] + & (x[:, target_q] ^ z[:, target_q]) + & (z[:, control_q] ^ ~x[:, target_q]) + ) + ^ ((x[:, target_q] ^ x[:, control_q]) & (z[:, target_q] ^ x[:, target_q])) + ) + x_target_q = x[:, control_q] ^ x[:, target_q] + z_control_q = z[:, control_q] ^ z[:, target_q] ^ x[:, target_q] + z_target_q = z[:, target_q] ^ x[:, control_q] + symplectic_matrix[:, target_q] = x_target_q + symplectic_matrix[:, nqubits + control_q] = z_control_q + symplectic_matrix[:, nqubits + target_q] = z_target_q + return symplectic_matrix + + +def CRX(symplectic_matrix, control_q, target_q, nqubits, theta): + # theta = 4 * n * pi + if theta % (4 * np.pi) == 0: + return I(symplectic_matrix, target_q, nqubits) + # theta = pi + 4 * n * pi + elif (theta / np.pi - 1) % 4 == 0: + symplectic_matrix = X(symplectic_matrix, target_q, nqubits) + symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) + symplectic_matrix = X(symplectic_matrix, target_q, nqubits) + return CY(symplectic_matrix, control_q, target_q, nqubits) + # theta = 2 * pi + 4 * n * pi + elif (theta / (2 * np.pi) - 1) % 2 == 0: + symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) + symplectic_matrix = Y(symplectic_matrix, target_q, nqubits) + symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) + return Y(symplectic_matrix, target_q, nqubits) + # theta = 3 * pi + 4 * n * pi + elif (theta / np.pi - 3) % 4 == 0: + symplectic_matrix = X(symplectic_matrix, target_q, nqubits) + symplectic_matrix = CY(symplectic_matrix, control_q, target_q, nqubits) + symplectic_matrix = X(symplectic_matrix, target_q, nqubits) + return CZ(symplectic_matrix, control_q, target_q, nqubits) + + +def CRZ(symplectic_matrix, control_q, target_q, nqubits, theta): + # theta = 4 * n * pi + if theta % (4 * np.pi) == 0: + return I(symplectic_matrix, target_q, nqubits) + # theta = pi + 4 * n * pi + elif (theta / np.pi - 1) % 4 == 0: + symplectic_matrix = X(symplectic_matrix, target_q, nqubits) + symplectic_matrix = CY(symplectic_matrix, control_q, target_q, nqubits) + symplectic_matrix = X(symplectic_matrix, target_q, nqubits) + return CNOT(symplectic_matrix, control_q, target_q, nqubits) + # theta = 2 * pi + 4 * n * pi + elif (theta / (2 * np.pi) - 1) % 2 == 0: + symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) + symplectic_matrix = X(symplectic_matrix, target_q, nqubits) + symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) + return X(symplectic_matrix, target_q, nqubits) + # theta = 3 * pi + 4 * n * pi + elif (theta / np.pi - 3) % 4 == 0: + symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) + symplectic_matrix = X(symplectic_matrix, target_q, nqubits) + symplectic_matrix = CY(symplectic_matrix, control_q, target_q, nqubits) + return X(symplectic_matrix, target_q, nqubits) + + +def CRY(symplectic_matrix, control_q, target_q, nqubits, theta): + # theta = 4 * n * pi + if theta % (4 * np.pi) == 0: + return I(symplectic_matrix, target_q, nqubits) + # theta = pi + 4 * n * pi + elif (theta / np.pi - 1) % 4 == 0: + symplectic_matrix = Z(symplectic_matrix, target_q, nqubits) + symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) + symplectic_matrix = Z(symplectic_matrix, target_q, nqubits) + return CZ(symplectic_matrix, control_q, target_q, nqubits) + # theta = 2 * pi + 4 * n * pi + elif (theta / (2 * np.pi) - 1) % 2 == 0: + return CRZ(symplectic_matrix, control_q, target_q, nqubits, theta) + # theta = 3 * pi + 4 * n * pi + elif (theta / np.pi - 3) % 4 == 0: + symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) + symplectic_matrix = Z(symplectic_matrix, target_q, nqubits) + symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) + return Z(symplectic_matrix, target_q, nqubits) + + +def ECR(symplectic_matrix, control_q, target_q, nqubits): + symplectic_matrix = S(symplectic_matrix, control_q, nqubits) + symplectic_matrix = SX(symplectic_matrix, target_q, nqubits) + symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) + return X(symplectic_matrix, control_q, nqubits) + + +def _exponent( + x1: np.ndarray, z1: np.ndarray, x2: np.ndarray, z2: np.ndarray +) -> np.ndarray: + """Helper function that computes the exponent to which i is raised for the product of the x and z paulis encoded in the symplectic matrix. This is used in _rowsum. The computation is performed parallely over the separated paulis x1[i], z1[i], x2[i] and z2[i]. + + Args: + x1 (np.array): Bits of the first x paulis. + z1 (np.array): Bits of the first z paulis. + x2 (np.array): Bits of the second x paulis. + z2 (np.array): Bits of the second z paulis. + + Returns: + ndarray: The calculated exponents. + """ + # this cannot be performed in the packed representation for measurements (thus packed rows) + # because bitwise arithmetic difference and sum are needed, which cannot be done directly + # in the packed representation. + return 2 * (x1 * x2 * (z2 - z1) + z1 * z2 * (x1 - x2)) - x1 * z2 + x2 * z1 + + +def _rowsum(symplectic_matrix, h, i, nqubits, determined=False): + """Helper function that updates the symplectic matrix by setting the h-th generator equal to the (i+h)-th one. This is done to keep track of the phase of the h-th row of the symplectic matrix (r[h]). The function is applied parallely over all the rows h and i passed. + + Args: + symplectic_matrix (np.array): Input symplectic matrix. + h (np.array): Indices of the rows encoding the generators to update. + i (np.array): Indices of the rows encoding the generators to use. + nqubits (int): Total number of qubits. + + Returns: + ndarray: The updated symplectic matrix. + """ + # calculate the exponent in the unpacked representation + xi, zi = symplectic_matrix[i, :nqubits], symplectic_matrix[i, nqubits:-1] + xh, zh = symplectic_matrix[h, :nqubits], symplectic_matrix[h, nqubits:-1] + exponents = _exponent(xi, zi, xh, zh) + ind = ( + 2 * symplectic_matrix[h, -1] + + 2 * symplectic_matrix[i, -1] + + np.sum(exponents, axis=-1) + ) % 4 == 0 + r = np.ones(h.shape[0], dtype=np.uint8) + r[ind] = 0 + + # the rest can be done in the packed representation + symplectic_matrix = _pack_for_measurements(symplectic_matrix, nqubits) + packed_n = _packed_size(nqubits) + xi, zi = symplectic_matrix[i, :packed_n], symplectic_matrix[i, packed_n:-1] + xh, zh = symplectic_matrix[h, :packed_n], symplectic_matrix[h, packed_n:-1] + xi_xh = xi ^ xh + zi_zh = zi ^ zh + if determined: + r = reduce(np.logical_xor, r) + xi_xh = reduce(np.logical_xor, xi_xh) + zi_zh = reduce(np.logical_xor, zi_zh) + symplectic_matrix[h[0], -1] = r + symplectic_matrix[h[0], :packed_n] = xi_xh + symplectic_matrix[h[0], packed_n:-1] = zi_zh + else: + symplectic_matrix[h, -1] = r + symplectic_matrix[h, :packed_n] = xi_xh + symplectic_matrix[h, packed_n:-1] = zi_zh + return _unpack_for_measurements(symplectic_matrix, nqubits) + + +def _determined_outcome(state, q, nqubits): + """Extracts the outcome for a measurement in case it is determined.""" + state[-1, :] = 0 + idx = (state[:nqubits, q].nonzero()[0] + nqubits).astype(np.uint) + if len(idx) == 0: + return state, state[-1, -1] + state = _rowsum( + state, + _dim_xz(nqubits) * np.ones(idx.shape, dtype=np.uint), + idx, + nqubits, + True, + ) + return state, state[-1, -1] + + +def _random_outcome(state, p, q, nqubits, outcome): + """Extracts the outcome for a measurement in case it is random.""" + p = p[0] + nqubits + state[p, q] = 0 + h = state[:-1, q].nonzero()[0] + state[p, q] = 1 + if h.shape[0] > 0: + state = _rowsum( + state, + h.astype(np.uint), + np.uint(p) * np.ones(h.shape[0], dtype=np.uint), + nqubits, + False, + ) + state[p - nqubits, :] = state[p, :] + #outcome = np.random.randint(2, size=1).item() + state[p, :] = 0 + state[p, -1] = outcome + state[p, nqubits + q] = 1 + return state, outcome + + +@cache +def _dim(nqubits): + """Returns the dimension of the symplectic matrix for a given number of qubits.""" + return _dim_xz(nqubits) + 1 + + +@cache +def _dim_xz(nqubits): + """Returns the dimension of the symplectic matrix (only the de/stabilizers generators part, + without the phases and scratch row) for a given number of qubits.""" + return 2 * nqubits + + +@cache +def _packed_size(n): + """Returns the size of an array of `n` booleans after packing.""" + return np.ceil(n / 8).astype(int) + + +def _packbits(array, axis): + return np.packbits(array, axis=axis) + + +def _unpackbits(array, axis, count): + return np.unpackbits(array, axis=axis, count=count) + + +def _pack_for_measurements(state, nqubits): + """Prepares the state for measurements by packing the rows of the X and Z sections of the symplectic matrix.""" + r, x, z = _get_rxz(state, nqubits) + x = _packbits(x, axis=1) + z = _packbits(z, axis=1) + return np.hstack((x, z, r[:, None])) + + +def _unpack_for_measurements(state, nqubits): + """Unpacks the symplectc matrix that was packed for measurements.""" + x = _unpackbits(state[:, : _packed_size(nqubits)], axis=1, count=nqubits) + z = _unpackbits(state[:, _packed_size(nqubits) : -1], axis=1, count=nqubits) + return np.hstack((x, z, state[:, -1][:, None])) + + +def _init_state_for_measurements(state, nqubits, collapse): + if collapse: + return _unpackbits(state, axis=0, count=_dim(nqubits)) + return state.copy() + +def _sample_random_outcomes(n_samples): + return np.random.randint(2, size=n_samples).tolist() + +def _M(state, qubits, nqubits, collapse=False): + sample = [] + state = _init_state_for_measurements(state, nqubits, collapse) + possible_outcomes = [] + for q in qubits: + p = state[nqubits:-1, q].nonzero()[0] + # random outcome, affects the state + if len(p) > 0: + state, _ = _random_outcome(state, p, q, nqubits, 0) + outcome = None + # determined outcome, state unchanged + else: + _, outcome = _determined_outcome(state, q, nqubits) + outcome = int(outcome) + sample.append(outcome) + if collapse: + state = _packbits(state, axis=0) + return sample + + +# valid for a standard basis measurement only +def M(state, qubits, nqubits, collapse=False, nshots=1): + if collapse and nshots != 1: + raise_error( + RuntimeError, + "Cannot generate multiple shots with a collapsing measurement.", + ) + sample = _M(state, qubits, nqubits, collapse) + samples = [ + ( + _sample_random_outcomes(nshots) + if outcome is None + else nshots * [outcome] + ) + for outcome in sample + ] + samples = list(zip(*samples)) + return samples + + +def cast(x, dtype=None, copy: bool = False): + if dtype is None: + dtype = "complex128" + if isinstance(x, np.ndarray): + return x.astype(dtype, copy=copy) + elif sparse.issparse(x): # pragma: no cover + return x.astype(dtype, copy=copy) + return np.asarray(x, dtype=dtype, copy=copy if copy else None) + + +def _clifford_pre_execution_reshape(state): + """Reshape and packing applied to the symplectic matrix before execution to prepare the state in the form needed by each engine. + + Args: + state (np.array): Input state. + + Returns: + (np.array) The packed and reshaped state. + """ + return _packbits(state, axis=0) + + +def _clifford_post_execution_reshape(state, nqubits: int): + """Reshape and unpacking applied to the state after execution to retrieve the standard symplectic matrix form. + + Args: + state (np.array): Input state. + nqubits (int): Number of qubits. + + Returns: + (np.array) The unpacked and reshaped state. + """ + state = _unpackbits(state, axis=0, count=_dim(nqubits))[: _dim(nqubits)] + return state + + +def identity_density_matrix(nqubits, normalize: bool = True): + state = np.eye(2**nqubits, dtype="complex128") + if normalize is True: # pragma: no cover + state /= 2**nqubits + return state + + +def set_seed(seed): + np.random.seed(seed) diff --git a/src/qibo/backends/.#_clifford_operations.py b/src/qibo/backends/.#_clifford_operations.py new file mode 120000 index 0000000000..c543502c6a --- /dev/null +++ b/src/qibo/backends/.#_clifford_operations.py @@ -0,0 +1 @@ +andrea@ubuntu-desktop.20019:1764055565 \ No newline at end of file diff --git a/src/qibo/backends/_clifford_operations.py b/src/qibo/backends/_clifford_operations.py index 3f79db813d..19a66878f8 100644 --- a/src/qibo/backends/_clifford_operations.py +++ b/src/qibo/backends/_clifford_operations.py @@ -426,7 +426,7 @@ def _determined_outcome(state, q, nqubits): return state, state[-1, -1] -def _random_outcome(state, p, q, nqubits): +def _random_outcome(state, p, q, nqubits, outcome): """Extracts the outcome for a measurement in case it is random.""" p = p[0] + nqubits state[p, q] = 0 @@ -441,7 +441,7 @@ def _random_outcome(state, p, q, nqubits): False, ) state[p - nqubits, :] = state[p, :] - outcome = np.random.randint(2, size=1).item() + # outcome = np.random.randint(2, size=1).item() state[p, :] = 0 state[p, -1] = outcome state[p, nqubits + q] = 1 @@ -496,10 +496,14 @@ def _init_state_for_measurements(state, nqubits, collapse): return state.copy() +def _sample_random_outcomes(n_samples): + return np.random.randint(2, size=n_samples).tolist() + + def _M(state, qubits, nqubits, collapse=False): sample = [] state = _init_state_for_measurements(state, nqubits, collapse) - # TODO: parallelize this and get rid of the loop + possible_outcomes = [] for q in qubits: p = state[nqubits:-1, q].nonzero()[0] # random outcome, affects the state @@ -525,11 +529,7 @@ def M(state, qubits, nqubits, collapse=False, nshots=1): ) sample = _M(state, qubits, nqubits, collapse) samples = [ - ( - np.random.randint(2, size=nshots).tolist() - if outcome is None - else nshots * [outcome] - ) + (_sample_random_outcomes(nshots) if outcome is None else nshots * [outcome]) for outcome in sample ] samples = list(zip(*samples)) diff --git a/tests/test_backends_clifford.py b/tests/test_backends_clifford.py index d333849f78..d62be0465e 100644 --- a/tests/test_backends_clifford.py +++ b/tests/test_backends_clifford.py @@ -161,7 +161,8 @@ def test_two_qubits_gates(backend, gate): [ 15, 25, - ], + ] + + list(range(100)), ) def test_random_clifford_circuit(backend, prob_qubits, binary, seed): np.random.seed(seed) From e809e86177560afc6216e8ff23a8321bb68bcc67 Mon Sep 17 00:00:00 2001 From: BrunoLiegiBastonLiegi Date: Tue, 25 Nov 2025 12:16:47 +0100 Subject: [PATCH 3/4] not working actually --- src/qibo/backends/#_clifford_operations.py# | 585 -------------------- src/qibo/backends/.#_clifford_operations.py | 1 - src/qibo/backends/_clifford_operations.py | 12 +- 3 files changed, 7 insertions(+), 591 deletions(-) delete mode 100644 src/qibo/backends/#_clifford_operations.py# delete mode 120000 src/qibo/backends/.#_clifford_operations.py diff --git a/src/qibo/backends/#_clifford_operations.py# b/src/qibo/backends/#_clifford_operations.py# deleted file mode 100644 index 786c2e7e4e..0000000000 --- a/src/qibo/backends/#_clifford_operations.py# +++ /dev/null @@ -1,585 +0,0 @@ -from functools import cache, reduce - -import numpy as np -from scipy import sparse - -from qibo.config import raise_error - -name = "numpy" - - -def _get_rxz(symplectic_matrix, nqubits): - return ( - symplectic_matrix[:, -1], - symplectic_matrix[:, :nqubits], - symplectic_matrix[:, nqubits:-1], - ) - - -def I(symplectic_matrix, q, nqubits): - return symplectic_matrix - - -def H(symplectic_matrix, q, nqubits): - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = r ^ (x[:, q] & z[:, q]) - symplectic_matrix[:, [q, nqubits + q]] = symplectic_matrix[:, [nqubits + q, q]] - return symplectic_matrix - - -def CNOT(symplectic_matrix, control_q, target_q, nqubits): - ind_zt = nqubits + target_q - ind_zc = nqubits + control_q - r = symplectic_matrix[:, -1] - xcq = symplectic_matrix[:, control_q] - xtq = symplectic_matrix[:, target_q] - ztq = symplectic_matrix[:, ind_zt] - zcq = symplectic_matrix[:, ind_zc] - symplectic_matrix[:, -1] = r ^ (xcq & ztq) & (xtq ^ ~zcq) - symplectic_matrix[:, target_q] = xtq ^ xcq - symplectic_matrix[:, ind_zc] = zcq ^ ztq - return symplectic_matrix - - -def CZ(symplectic_matrix, control_q, target_q, nqubits): - """Decomposition --> H-CNOT-H""" - ind_zt = nqubits + target_q - ind_zc = nqubits + control_q - r = symplectic_matrix[:, -1] - xcq = symplectic_matrix[:, control_q] - xtq = symplectic_matrix[:, target_q] - ztq = symplectic_matrix[:, ind_zt] - zcq = symplectic_matrix[:, ind_zc] - ztq_xor_xcq = ztq ^ xcq - symplectic_matrix[:, -1] = ( - r ^ (xtq & ztq) ^ (xcq & xtq & (ztq ^ ~zcq)) ^ (xtq & ztq_xor_xcq) - ) - z_control_q = xtq ^ zcq - z_target_q = ztq_xor_xcq - symplectic_matrix[:, ind_zc] = z_control_q - symplectic_matrix[:, ind_zt] = z_target_q - return symplectic_matrix - - -def S(symplectic_matrix, q, nqubits): - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = r ^ (x[:, q] & z[:, q]) - symplectic_matrix[:, q + nqubits] = z[:, q] ^ x[:, q] - return symplectic_matrix - - -def Z(symplectic_matrix, q, nqubits): - """Decomposition --> S-S""" - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = r ^ ((x[:, q] & z[:, q]) ^ x[:, q] & (z[:, q] ^ x[:, q])) - return symplectic_matrix - - -def X(symplectic_matrix, q, nqubits): - """Decomposition --> H-S-S-H""" - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = r ^ (z[:, q] & (z[:, q] ^ x[:, q])) ^ (z[:, q] & x[:, q]) - return symplectic_matrix - - -def Y(symplectic_matrix, q, nqubits): - """Decomposition --> S-S-H-S-S-H""" - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = ( - r ^ (z[:, q] & (z[:, q] ^ x[:, q])) ^ (x[:, q] & (z[:, q] ^ x[:, q])) - ) - return symplectic_matrix - - -def SX(symplectic_matrix, q, nqubits): - """Decomposition --> H-S-H""" - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = r ^ (z[:, q] & (z[:, q] ^ x[:, q])) - symplectic_matrix[:, q] = z[:, q] ^ x[:, q] - return symplectic_matrix - - -def SDG(symplectic_matrix, q, nqubits): - """Decomposition --> S-S-S""" - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = r ^ (x[:, q] & (z[:, q] ^ x[:, q])) - symplectic_matrix[:, nqubits + q] = z[:, q] ^ x[:, q] - return symplectic_matrix - - -def SXDG(symplectic_matrix, q, nqubits): - """Decomposition --> H-S-S-S-H""" - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = r ^ (z[:, q] & x[:, q]) - symplectic_matrix[:, q] = z[:, q] ^ x[:, q] - return symplectic_matrix - - -def RX(symplectic_matrix, q, nqubits, theta): - if theta % (2 * np.pi) == 0: - return I(symplectic_matrix, q, nqubits) - elif (theta / np.pi - 1) % 2 == 0: - return X(symplectic_matrix, q, nqubits) - elif (theta / (np.pi / 2) - 1) % 4 == 0: - return SX(symplectic_matrix, q, nqubits) - else: # theta == 3*pi/2 + 2*n*pi - return SXDG(symplectic_matrix, q, nqubits) - - -def RZ(symplectic_matrix, q, nqubits, theta): - if theta % (2 * np.pi) == 0: - return I(symplectic_matrix, q, nqubits) - elif (theta / np.pi - 1) % 2 == 0: - return Z(symplectic_matrix, q, nqubits) - elif (theta / (np.pi / 2) - 1) % 4 == 0: - return S(symplectic_matrix, q, nqubits) - else: # theta == 3*pi/2 + 2*n*pi - return SDG(symplectic_matrix, q, nqubits) - - -def RY_pi(symplectic_matrix, q, nqubits): - """Decomposition --> H-S-S""" - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = r ^ (x[:, q] & (z[:, q] ^ x[:, q])) - symplectic_matrix[:, [nqubits + q, q]] = symplectic_matrix[:, [q, nqubits + q]] - return symplectic_matrix - - -def RY_3pi_2(symplectic_matrix, q, nqubits): - """Decomposition --> H-S-S-H-S-S-H-S-S""" - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = r ^ (z[:, q] & (z[:, q] ^ x[:, q])) - symplectic_matrix[:, [nqubits + q, q]] = symplectic_matrix[:, [q, nqubits + q]] - return symplectic_matrix - - -def RY(symplectic_matrix, q, nqubits, theta): - if theta % (2 * np.pi) == 0: - return I(symplectic_matrix, q, nqubits) - elif (theta / np.pi - 1) % 2 == 0: - return Y(symplectic_matrix, q, nqubits) - elif (theta / (np.pi / 2) - 1) % 4 == 0: - """Decomposition --> H-S-S""" - return RY_pi(symplectic_matrix, q, nqubits) - else: # theta == 3*pi/2 + 2*n*pi - """Decomposition --> H-S-S-H-S-S-H-S-S""" - return RY_3pi_2(symplectic_matrix, q, nqubits) - - -def GPI2(symplectic_matrix, q, nqubits, phi): - if phi % (2 * np.pi) == 0: - return RX(symplectic_matrix, q, nqubits, np.pi / 2) - - if (phi / np.pi - 1) % 2 == 0: - return RX(symplectic_matrix, q, nqubits, -np.pi / 2) - - if (phi / (np.pi / 2) - 1) % 4 == 0: - return RY(symplectic_matrix, q, nqubits, np.pi / 2) - - # theta == 3*pi/2 + 2*n*pi - return RY(symplectic_matrix, q, nqubits, -np.pi / 2) - - -def SWAP(symplectic_matrix, control_q, target_q, nqubits): - """Decomposition --> CNOT-CNOT-CNOT""" - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = ( - r - ^ (x[:, control_q] & z[:, target_q] & (x[:, target_q] ^ ~z[:, control_q])) - ^ ( - (x[:, target_q] ^ x[:, control_q]) - & (z[:, target_q] ^ z[:, control_q]) - & (z[:, target_q] ^ ~x[:, control_q]) - ) - ^ ( - x[:, target_q] - & z[:, control_q] - & (x[:, control_q] ^ x[:, target_q] ^ z[:, control_q] ^ ~z[:, target_q]) - ) - ) - symplectic_matrix[ - :, [control_q, target_q, nqubits + control_q, nqubits + target_q] - ] = symplectic_matrix[ - :, [target_q, control_q, nqubits + target_q, nqubits + control_q] - ] - return symplectic_matrix - - -def iSWAP(symplectic_matrix, control_q, target_q, nqubits): - """Decomposition --> H-CNOT-CNOT-H-S-S""" - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = ( - r - ^ (x[:, target_q] & z[:, target_q]) - ^ (x[:, control_q] & z[:, control_q]) - ^ (x[:, control_q] & (z[:, control_q] ^ x[:, control_q])) - ^ ( - (z[:, control_q] ^ x[:, control_q]) - & (z[:, target_q] ^ x[:, target_q]) - & (x[:, target_q] ^ ~x[:, control_q]) - ) - ^ ( - (x[:, target_q] ^ z[:, control_q] ^ x[:, control_q]) - & (x[:, target_q] ^ z[:, target_q] ^ x[:, control_q]) - & (x[:, target_q] ^ z[:, target_q] ^ x[:, control_q] ^ ~z[:, control_q]) - ) - ^ (x[:, control_q] & (x[:, target_q] ^ x[:, control_q] ^ z[:, control_q])) - ) - z_control_q = x[:, target_q] ^ z[:, target_q] ^ x[:, control_q] - z_target_q = x[:, target_q] ^ z[:, control_q] ^ x[:, control_q] - symplectic_matrix[:, nqubits + control_q] = z_control_q - symplectic_matrix[:, nqubits + target_q] = z_target_q - symplectic_matrix[:, [control_q, target_q]] = symplectic_matrix[ - :, [target_q, control_q] - ] - return symplectic_matrix - - -def FSWAP(symplectic_matrix, control_q, target_q, nqubits): - """Decomposition --> X-CNOT-RY-CNOT-RY-CNOT-CNOT-X""" - symplectic_matrix = X(symplectic_matrix, target_q, nqubits) - symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) - symplectic_matrix = RY(symplectic_matrix, control_q, nqubits, np.pi / 2) - symplectic_matrix = CNOT(symplectic_matrix, target_q, control_q, nqubits) - symplectic_matrix = RY(symplectic_matrix, control_q, nqubits, -np.pi / 2) - symplectic_matrix = CNOT(symplectic_matrix, target_q, control_q, nqubits) - symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) - return X(symplectic_matrix, control_q, nqubits) - - -def CY(symplectic_matrix, control_q, target_q, nqubits): - """Decomposition --> S-CNOT-SDG""" - r, x, z = _get_rxz(symplectic_matrix, nqubits) - symplectic_matrix[:, -1] = ( - r - ^ (x[:, target_q] & (z[:, target_q] ^ x[:, target_q])) - ^ ( - x[:, control_q] - & (x[:, target_q] ^ z[:, target_q]) - & (z[:, control_q] ^ ~x[:, target_q]) - ) - ^ ((x[:, target_q] ^ x[:, control_q]) & (z[:, target_q] ^ x[:, target_q])) - ) - x_target_q = x[:, control_q] ^ x[:, target_q] - z_control_q = z[:, control_q] ^ z[:, target_q] ^ x[:, target_q] - z_target_q = z[:, target_q] ^ x[:, control_q] - symplectic_matrix[:, target_q] = x_target_q - symplectic_matrix[:, nqubits + control_q] = z_control_q - symplectic_matrix[:, nqubits + target_q] = z_target_q - return symplectic_matrix - - -def CRX(symplectic_matrix, control_q, target_q, nqubits, theta): - # theta = 4 * n * pi - if theta % (4 * np.pi) == 0: - return I(symplectic_matrix, target_q, nqubits) - # theta = pi + 4 * n * pi - elif (theta / np.pi - 1) % 4 == 0: - symplectic_matrix = X(symplectic_matrix, target_q, nqubits) - symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) - symplectic_matrix = X(symplectic_matrix, target_q, nqubits) - return CY(symplectic_matrix, control_q, target_q, nqubits) - # theta = 2 * pi + 4 * n * pi - elif (theta / (2 * np.pi) - 1) % 2 == 0: - symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) - symplectic_matrix = Y(symplectic_matrix, target_q, nqubits) - symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) - return Y(symplectic_matrix, target_q, nqubits) - # theta = 3 * pi + 4 * n * pi - elif (theta / np.pi - 3) % 4 == 0: - symplectic_matrix = X(symplectic_matrix, target_q, nqubits) - symplectic_matrix = CY(symplectic_matrix, control_q, target_q, nqubits) - symplectic_matrix = X(symplectic_matrix, target_q, nqubits) - return CZ(symplectic_matrix, control_q, target_q, nqubits) - - -def CRZ(symplectic_matrix, control_q, target_q, nqubits, theta): - # theta = 4 * n * pi - if theta % (4 * np.pi) == 0: - return I(symplectic_matrix, target_q, nqubits) - # theta = pi + 4 * n * pi - elif (theta / np.pi - 1) % 4 == 0: - symplectic_matrix = X(symplectic_matrix, target_q, nqubits) - symplectic_matrix = CY(symplectic_matrix, control_q, target_q, nqubits) - symplectic_matrix = X(symplectic_matrix, target_q, nqubits) - return CNOT(symplectic_matrix, control_q, target_q, nqubits) - # theta = 2 * pi + 4 * n * pi - elif (theta / (2 * np.pi) - 1) % 2 == 0: - symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) - symplectic_matrix = X(symplectic_matrix, target_q, nqubits) - symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) - return X(symplectic_matrix, target_q, nqubits) - # theta = 3 * pi + 4 * n * pi - elif (theta / np.pi - 3) % 4 == 0: - symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) - symplectic_matrix = X(symplectic_matrix, target_q, nqubits) - symplectic_matrix = CY(symplectic_matrix, control_q, target_q, nqubits) - return X(symplectic_matrix, target_q, nqubits) - - -def CRY(symplectic_matrix, control_q, target_q, nqubits, theta): - # theta = 4 * n * pi - if theta % (4 * np.pi) == 0: - return I(symplectic_matrix, target_q, nqubits) - # theta = pi + 4 * n * pi - elif (theta / np.pi - 1) % 4 == 0: - symplectic_matrix = Z(symplectic_matrix, target_q, nqubits) - symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) - symplectic_matrix = Z(symplectic_matrix, target_q, nqubits) - return CZ(symplectic_matrix, control_q, target_q, nqubits) - # theta = 2 * pi + 4 * n * pi - elif (theta / (2 * np.pi) - 1) % 2 == 0: - return CRZ(symplectic_matrix, control_q, target_q, nqubits, theta) - # theta = 3 * pi + 4 * n * pi - elif (theta / np.pi - 3) % 4 == 0: - symplectic_matrix = CZ(symplectic_matrix, control_q, target_q, nqubits) - symplectic_matrix = Z(symplectic_matrix, target_q, nqubits) - symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) - return Z(symplectic_matrix, target_q, nqubits) - - -def ECR(symplectic_matrix, control_q, target_q, nqubits): - symplectic_matrix = S(symplectic_matrix, control_q, nqubits) - symplectic_matrix = SX(symplectic_matrix, target_q, nqubits) - symplectic_matrix = CNOT(symplectic_matrix, control_q, target_q, nqubits) - return X(symplectic_matrix, control_q, nqubits) - - -def _exponent( - x1: np.ndarray, z1: np.ndarray, x2: np.ndarray, z2: np.ndarray -) -> np.ndarray: - """Helper function that computes the exponent to which i is raised for the product of the x and z paulis encoded in the symplectic matrix. This is used in _rowsum. The computation is performed parallely over the separated paulis x1[i], z1[i], x2[i] and z2[i]. - - Args: - x1 (np.array): Bits of the first x paulis. - z1 (np.array): Bits of the first z paulis. - x2 (np.array): Bits of the second x paulis. - z2 (np.array): Bits of the second z paulis. - - Returns: - ndarray: The calculated exponents. - """ - # this cannot be performed in the packed representation for measurements (thus packed rows) - # because bitwise arithmetic difference and sum are needed, which cannot be done directly - # in the packed representation. - return 2 * (x1 * x2 * (z2 - z1) + z1 * z2 * (x1 - x2)) - x1 * z2 + x2 * z1 - - -def _rowsum(symplectic_matrix, h, i, nqubits, determined=False): - """Helper function that updates the symplectic matrix by setting the h-th generator equal to the (i+h)-th one. This is done to keep track of the phase of the h-th row of the symplectic matrix (r[h]). The function is applied parallely over all the rows h and i passed. - - Args: - symplectic_matrix (np.array): Input symplectic matrix. - h (np.array): Indices of the rows encoding the generators to update. - i (np.array): Indices of the rows encoding the generators to use. - nqubits (int): Total number of qubits. - - Returns: - ndarray: The updated symplectic matrix. - """ - # calculate the exponent in the unpacked representation - xi, zi = symplectic_matrix[i, :nqubits], symplectic_matrix[i, nqubits:-1] - xh, zh = symplectic_matrix[h, :nqubits], symplectic_matrix[h, nqubits:-1] - exponents = _exponent(xi, zi, xh, zh) - ind = ( - 2 * symplectic_matrix[h, -1] - + 2 * symplectic_matrix[i, -1] - + np.sum(exponents, axis=-1) - ) % 4 == 0 - r = np.ones(h.shape[0], dtype=np.uint8) - r[ind] = 0 - - # the rest can be done in the packed representation - symplectic_matrix = _pack_for_measurements(symplectic_matrix, nqubits) - packed_n = _packed_size(nqubits) - xi, zi = symplectic_matrix[i, :packed_n], symplectic_matrix[i, packed_n:-1] - xh, zh = symplectic_matrix[h, :packed_n], symplectic_matrix[h, packed_n:-1] - xi_xh = xi ^ xh - zi_zh = zi ^ zh - if determined: - r = reduce(np.logical_xor, r) - xi_xh = reduce(np.logical_xor, xi_xh) - zi_zh = reduce(np.logical_xor, zi_zh) - symplectic_matrix[h[0], -1] = r - symplectic_matrix[h[0], :packed_n] = xi_xh - symplectic_matrix[h[0], packed_n:-1] = zi_zh - else: - symplectic_matrix[h, -1] = r - symplectic_matrix[h, :packed_n] = xi_xh - symplectic_matrix[h, packed_n:-1] = zi_zh - return _unpack_for_measurements(symplectic_matrix, nqubits) - - -def _determined_outcome(state, q, nqubits): - """Extracts the outcome for a measurement in case it is determined.""" - state[-1, :] = 0 - idx = (state[:nqubits, q].nonzero()[0] + nqubits).astype(np.uint) - if len(idx) == 0: - return state, state[-1, -1] - state = _rowsum( - state, - _dim_xz(nqubits) * np.ones(idx.shape, dtype=np.uint), - idx, - nqubits, - True, - ) - return state, state[-1, -1] - - -def _random_outcome(state, p, q, nqubits, outcome): - """Extracts the outcome for a measurement in case it is random.""" - p = p[0] + nqubits - state[p, q] = 0 - h = state[:-1, q].nonzero()[0] - state[p, q] = 1 - if h.shape[0] > 0: - state = _rowsum( - state, - h.astype(np.uint), - np.uint(p) * np.ones(h.shape[0], dtype=np.uint), - nqubits, - False, - ) - state[p - nqubits, :] = state[p, :] - #outcome = np.random.randint(2, size=1).item() - state[p, :] = 0 - state[p, -1] = outcome - state[p, nqubits + q] = 1 - return state, outcome - - -@cache -def _dim(nqubits): - """Returns the dimension of the symplectic matrix for a given number of qubits.""" - return _dim_xz(nqubits) + 1 - - -@cache -def _dim_xz(nqubits): - """Returns the dimension of the symplectic matrix (only the de/stabilizers generators part, - without the phases and scratch row) for a given number of qubits.""" - return 2 * nqubits - - -@cache -def _packed_size(n): - """Returns the size of an array of `n` booleans after packing.""" - return np.ceil(n / 8).astype(int) - - -def _packbits(array, axis): - return np.packbits(array, axis=axis) - - -def _unpackbits(array, axis, count): - return np.unpackbits(array, axis=axis, count=count) - - -def _pack_for_measurements(state, nqubits): - """Prepares the state for measurements by packing the rows of the X and Z sections of the symplectic matrix.""" - r, x, z = _get_rxz(state, nqubits) - x = _packbits(x, axis=1) - z = _packbits(z, axis=1) - return np.hstack((x, z, r[:, None])) - - -def _unpack_for_measurements(state, nqubits): - """Unpacks the symplectc matrix that was packed for measurements.""" - x = _unpackbits(state[:, : _packed_size(nqubits)], axis=1, count=nqubits) - z = _unpackbits(state[:, _packed_size(nqubits) : -1], axis=1, count=nqubits) - return np.hstack((x, z, state[:, -1][:, None])) - - -def _init_state_for_measurements(state, nqubits, collapse): - if collapse: - return _unpackbits(state, axis=0, count=_dim(nqubits)) - return state.copy() - -def _sample_random_outcomes(n_samples): - return np.random.randint(2, size=n_samples).tolist() - -def _M(state, qubits, nqubits, collapse=False): - sample = [] - state = _init_state_for_measurements(state, nqubits, collapse) - possible_outcomes = [] - for q in qubits: - p = state[nqubits:-1, q].nonzero()[0] - # random outcome, affects the state - if len(p) > 0: - state, _ = _random_outcome(state, p, q, nqubits, 0) - outcome = None - # determined outcome, state unchanged - else: - _, outcome = _determined_outcome(state, q, nqubits) - outcome = int(outcome) - sample.append(outcome) - if collapse: - state = _packbits(state, axis=0) - return sample - - -# valid for a standard basis measurement only -def M(state, qubits, nqubits, collapse=False, nshots=1): - if collapse and nshots != 1: - raise_error( - RuntimeError, - "Cannot generate multiple shots with a collapsing measurement.", - ) - sample = _M(state, qubits, nqubits, collapse) - samples = [ - ( - _sample_random_outcomes(nshots) - if outcome is None - else nshots * [outcome] - ) - for outcome in sample - ] - samples = list(zip(*samples)) - return samples - - -def cast(x, dtype=None, copy: bool = False): - if dtype is None: - dtype = "complex128" - if isinstance(x, np.ndarray): - return x.astype(dtype, copy=copy) - elif sparse.issparse(x): # pragma: no cover - return x.astype(dtype, copy=copy) - return np.asarray(x, dtype=dtype, copy=copy if copy else None) - - -def _clifford_pre_execution_reshape(state): - """Reshape and packing applied to the symplectic matrix before execution to prepare the state in the form needed by each engine. - - Args: - state (np.array): Input state. - - Returns: - (np.array) The packed and reshaped state. - """ - return _packbits(state, axis=0) - - -def _clifford_post_execution_reshape(state, nqubits: int): - """Reshape and unpacking applied to the state after execution to retrieve the standard symplectic matrix form. - - Args: - state (np.array): Input state. - nqubits (int): Number of qubits. - - Returns: - (np.array) The unpacked and reshaped state. - """ - state = _unpackbits(state, axis=0, count=_dim(nqubits))[: _dim(nqubits)] - return state - - -def identity_density_matrix(nqubits, normalize: bool = True): - state = np.eye(2**nqubits, dtype="complex128") - if normalize is True: # pragma: no cover - state /= 2**nqubits - return state - - -def set_seed(seed): - np.random.seed(seed) diff --git a/src/qibo/backends/.#_clifford_operations.py b/src/qibo/backends/.#_clifford_operations.py deleted file mode 120000 index c543502c6a..0000000000 --- a/src/qibo/backends/.#_clifford_operations.py +++ /dev/null @@ -1 +0,0 @@ -andrea@ubuntu-desktop.20019:1764055565 \ No newline at end of file diff --git a/src/qibo/backends/_clifford_operations.py b/src/qibo/backends/_clifford_operations.py index 19a66878f8..786c2e7e4e 100644 --- a/src/qibo/backends/_clifford_operations.py +++ b/src/qibo/backends/_clifford_operations.py @@ -441,7 +441,7 @@ def _random_outcome(state, p, q, nqubits, outcome): False, ) state[p - nqubits, :] = state[p, :] - # outcome = np.random.randint(2, size=1).item() + #outcome = np.random.randint(2, size=1).item() state[p, :] = 0 state[p, -1] = outcome state[p, nqubits + q] = 1 @@ -495,11 +495,9 @@ def _init_state_for_measurements(state, nqubits, collapse): return _unpackbits(state, axis=0, count=_dim(nqubits)) return state.copy() - def _sample_random_outcomes(n_samples): return np.random.randint(2, size=n_samples).tolist() - def _M(state, qubits, nqubits, collapse=False): sample = [] state = _init_state_for_measurements(state, nqubits, collapse) @@ -508,7 +506,7 @@ def _M(state, qubits, nqubits, collapse=False): p = state[nqubits:-1, q].nonzero()[0] # random outcome, affects the state if len(p) > 0: - state, _ = _random_outcome(state, p, q, nqubits) + state, _ = _random_outcome(state, p, q, nqubits, 0) outcome = None # determined outcome, state unchanged else: @@ -529,7 +527,11 @@ def M(state, qubits, nqubits, collapse=False, nshots=1): ) sample = _M(state, qubits, nqubits, collapse) samples = [ - (_sample_random_outcomes(nshots) if outcome is None else nshots * [outcome]) + ( + _sample_random_outcomes(nshots) + if outcome is None + else nshots * [outcome] + ) for outcome in sample ] samples = list(zip(*samples)) From a5ccb34f0012c7e87b9d86bdead102ac6958782f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 25 Nov 2025 11:18:19 +0000 Subject: [PATCH 4/4] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/qibo/backends/_clifford_operations.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/qibo/backends/_clifford_operations.py b/src/qibo/backends/_clifford_operations.py index 786c2e7e4e..efcb212852 100644 --- a/src/qibo/backends/_clifford_operations.py +++ b/src/qibo/backends/_clifford_operations.py @@ -441,7 +441,7 @@ def _random_outcome(state, p, q, nqubits, outcome): False, ) state[p - nqubits, :] = state[p, :] - #outcome = np.random.randint(2, size=1).item() + # outcome = np.random.randint(2, size=1).item() state[p, :] = 0 state[p, -1] = outcome state[p, nqubits + q] = 1 @@ -495,9 +495,11 @@ def _init_state_for_measurements(state, nqubits, collapse): return _unpackbits(state, axis=0, count=_dim(nqubits)) return state.copy() + def _sample_random_outcomes(n_samples): return np.random.randint(2, size=n_samples).tolist() + def _M(state, qubits, nqubits, collapse=False): sample = [] state = _init_state_for_measurements(state, nqubits, collapse) @@ -527,11 +529,7 @@ def M(state, qubits, nqubits, collapse=False, nshots=1): ) sample = _M(state, qubits, nqubits, collapse) samples = [ - ( - _sample_random_outcomes(nshots) - if outcome is None - else nshots * [outcome] - ) + (_sample_random_outcomes(nshots) if outcome is None else nshots * [outcome]) for outcome in sample ] samples = list(zip(*samples))