From 7717d7c57273224e9888d8ab7d8bd33e2f13ad30 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 6 Jul 2026 15:17:58 +0800 Subject: [PATCH 01/36] Add failed example to tests --- tests/test_measurement_optimisation.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/tests/test_measurement_optimisation.py b/tests/test_measurement_optimisation.py index c4f6a0cc..f007e695 100644 --- a/tests/test_measurement_optimisation.py +++ b/tests/test_measurement_optimisation.py @@ -19,6 +19,12 @@ def test_gc_measurement_mapping(): """Remaining coverage tests for _gc_measurement_mapping""" + # Check that qubits to measure remain unchanged + ham = SymbolicHamiltonian(X(2) * X(3) + Y(2) * Y(3) + Z(2) * Z(3), nqubits=4) + mapping, m_gates = _gc_measurement_mapping(ham.form, ham.nqubits, "chong") + assert {term.target_qubit for term in mapping.values() if hasattr(term, "target_qubit")} == {2, 3} + + # Single term Hamiltonian ham = SymbolicHamiltonian(Z(2)) mapping, m_gates = _gc_measurement_mapping(ham.form, 2, "izmaylov") assert mapping == {"Z2": ham.form} # Single term expression should remain unchanged From a3ad067ec7b2eab86876f88003df4a5d2e6e8760 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 6 Jul 2026 15:21:19 +0800 Subject: [PATCH 02/36] Add mapping back to original set of qubits --- src/qibochem/measurement/optimization.py | 15 +++++++++++++-- src/qibochem/measurement/util.py | 6 +++++- 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index cfa7acf5..4476a5a1 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -16,6 +16,7 @@ from qibochem.measurement.util import ( _binary_gaussian_elimination, _binary_nullspace, + _get_qubit, _get_sigma_terms, _group_commuting_terms, _lagrangian_subspace, @@ -124,6 +125,7 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl ] # Otherwise, expression is a sum of terms term_list = [_term_to_string(term) for term in expression.args if _term_to_string(term)[0] in ("X", "Y", "Z")] + term_qubits = sorted({_get_qubit(op) for term in term_list for op in term.split()}) v_subspace = np.array([_pauli_to_symplectic(terms.split(), nqubits) for terms in term_list]) v_basis = _binary_gaussian_elimination(v_subspace) @@ -135,24 +137,33 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl # Interchange the 1st/2nd half of the indices to get nullspace in a symplectic sense nullspace = np.concatenate((nullspace[:, dim_symplectic:], nullspace[:, :dim_symplectic]), axis=1) v_basis = _lagrangian_subspace(nullspace) + # Different methods of circuit synthesis if method == "chong": x_result = _solve_linear_system(v_basis, v_subspace) + # Map the solution onto the original set of qubits + qubit_map = { + q: initial_q for q, initial_q in zip(sorted({q for pauli_op in x_result for q in pauli_op}), term_qubits) + } phase_factors = [_phase_factor(v_basis[pauli_op]) for pauli_op in x_result] u_gates = _synthesise_circuit(v_basis) mapping = { - term: phase * prod(Z(_i) for _i in soln) for term, phase, soln in zip(term_list, phase_factors, x_result) + term: phase * prod(Z(qubit_map[_i]) for _i in soln) + for term, phase, soln in zip(term_list, phase_factors, x_result) } elif method == "izmaylov": v_basis = _sort_tau_terms(v_basis) new_tau_terms, sigma_terms = _get_sigma_terms(v_basis) x_result = _solve_linear_system(new_tau_terms, v_subspace) + qubit_map = { + q: initial_q for q, initial_q in zip(sorted({q for pauli_op in x_result for q in pauli_op}), term_qubits) + } phase_factors = [_phase_factor(new_tau_terms[pauli_op]) for pauli_op in x_result] tau_term_str = [_symplectic_to_pauli(tau_i) for tau_i in new_tau_terms] sigma_term_str = [_symplectic_to_pauli(sigma_i) for sigma_i in sigma_terms] qwc_terms = [_symplectic_to_pauli(sum(sigma_terms[_x] for _x in pauli_op)) for pauli_op in x_result] mapping = { - term: phase * prod([getattr(symbols, sigma[0])(int(sigma[1:])) for sigma in pauli_op]) + term: phase * prod([getattr(symbols, sigma[0])(qubit_map[int(sigma[1:])]) for sigma in pauli_op]) for term, phase, pauli_op in zip(term_list, phase_factors, qwc_terms) } # Define the measurement gates diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index b559c131..59d90cd3 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -185,7 +185,11 @@ def _lagrangian_subspace(vector_space: np.ndarray) -> np.ndarray: ) space_to_orthogonalize = space_to_orthogonalize % 2 - cp_vector_space = np.append([anticommuting_vectors[0]], space_to_orthogonalize, axis=0) + # Preferentially select Z over X + first_nonzero_col = np.argmax(anticommuting_vectors, axis=1) + selected_vector = anticommuting_vectors[np.argmax(first_nonzero_col)] + + cp_vector_space = np.append([selected_vector], space_to_orthogonalize, axis=0) cp_vector_space = _binary_gaussian_elimination(cp_vector_space) return cp_vector_space From 1fc0b224e0ae3f887cbcb40b0c5546cfb825eb3e Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 6 Jul 2026 15:39:00 +0800 Subject: [PATCH 03/36] Apply pylint suggestions --- src/qibochem/measurement/optimization.py | 8 ++------ tests/test_measurement_optimisation.py | 1 - 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index 4476a5a1..0267c607 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -142,9 +142,7 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl if method == "chong": x_result = _solve_linear_system(v_basis, v_subspace) # Map the solution onto the original set of qubits - qubit_map = { - q: initial_q for q, initial_q in zip(sorted({q for pauli_op in x_result for q in pauli_op}), term_qubits) - } + qubit_map = dict(zip(sorted({q for pauli_op in x_result for q in pauli_op}), term_qubits)) phase_factors = [_phase_factor(v_basis[pauli_op]) for pauli_op in x_result] u_gates = _synthesise_circuit(v_basis) mapping = { @@ -155,9 +153,7 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl v_basis = _sort_tau_terms(v_basis) new_tau_terms, sigma_terms = _get_sigma_terms(v_basis) x_result = _solve_linear_system(new_tau_terms, v_subspace) - qubit_map = { - q: initial_q for q, initial_q in zip(sorted({q for pauli_op in x_result for q in pauli_op}), term_qubits) - } + qubit_map = dict(zip(sorted({q for pauli_op in x_result for q in pauli_op}), term_qubits)) phase_factors = [_phase_factor(new_tau_terms[pauli_op]) for pauli_op in x_result] tau_term_str = [_symplectic_to_pauli(tau_i) for tau_i in new_tau_terms] sigma_term_str = [_symplectic_to_pauli(sigma_i) for sigma_i in sigma_terms] diff --git a/tests/test_measurement_optimisation.py b/tests/test_measurement_optimisation.py index f007e695..4e21060f 100644 --- a/tests/test_measurement_optimisation.py +++ b/tests/test_measurement_optimisation.py @@ -2,7 +2,6 @@ Test functionality to reduce the measurement cost of running VQE """ -import numpy as np import pytest from qibo.hamiltonians import SymbolicHamiltonian from qibo.symbols import X, Y, Z From b4dd22647895424a2a3e88c9ca4b0446c0bbf1c9 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Fri, 24 Jul 2026 16:10:30 +0800 Subject: [PATCH 04/36] Additional bug fix --- src/qibochem/measurement/optimization.py | 1 + src/qibochem/measurement/util.py | 7 ++++++- tests/test_expectation_samples.py | 19 +++++++++++++++++++ 3 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index 0267c607..2d3ade2b 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -136,6 +136,7 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl nullspace = _binary_nullspace(v_basis) # Interchange the 1st/2nd half of the indices to get nullspace in a symplectic sense nullspace = np.concatenate((nullspace[:, dim_symplectic:], nullspace[:, :dim_symplectic]), axis=1) + nullspace = _binary_gaussian_elimination(nullspace) v_basis = _lagrangian_subspace(nullspace) # Different methods of circuit synthesis diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index 59d90cd3..db2e8a72 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -175,7 +175,6 @@ def _lagrangian_subspace(vector_space: np.ndarray) -> np.ndarray: if cp_vector_space.shape[0] == (cp_vector_space.shape[1] // 2): break - # Remove the two anti-commuting vectors from the basis space_to_orthogonalize = np.delete(cp_vector_space, anticommuting_vector_indices, axis=0) for _i1, vector in enumerate(space_to_orthogonalize): @@ -202,6 +201,12 @@ def _sort_tau_terms(v_basis: np.ndarray) -> np.ndarray: will return [['X0', 'X2'], ['Z1'], ['Z0', 'Z2'], ['Z3', 'Z5'], ['Z4'], ['Z1', 'X3', 'Z4', 'X5']] """ + # Check to see if sorting needed + dim = v_basis.shape[0] + if all(v_basis[i, i] or v_basis[i, i + dim] for i in range(dim)): + return v_basis + + # TODO: Can refactor this? # Convert the basis set to strings for easier sorting pauli_terms = [_symplectic_to_pauli(vector) for vector in v_basis] dim = len(pauli_terms) diff --git a/tests/test_expectation_samples.py b/tests/test_expectation_samples.py index 4155f949..88808e66 100644 --- a/tests/test_expectation_samples.py +++ b/tests/test_expectation_samples.py @@ -2,6 +2,7 @@ Test expectation functionality """ +import numpy as np import pytest from qibo import Circuit, gates from qibo.hamiltonians import SymbolicHamiltonian @@ -112,6 +113,24 @@ def test_measurement_grouping_functionality(hamiltonian): assert test == pytest.approx(expected, abs=0.08) +def test_measurement_grouping_extra_tests(): + """Additional tests for some previously undiscovered bugs""" + for method in ("gc", "gc2"): + hamiltonian = SymbolicHamiltonian(0.5 * X(0) * Y(1) + Z(0) * Z(1) * Z(2), nqubits=3) + circuit = Circuit(3) + circuit.add(gates.H(0)) # = +1 + circuit.add(gates.RX(1, theta=-np.pi / 2)) # = +1 -> exact = +0.5 + result = expectation_from_samples(circuit, hamiltonian, n_shots=100_000, grouping=method) + assert result == pytest.approx(hamiltonian.expectation(circuit), abs=0.01) + + hamiltonian = SymbolicHamiltonian(0.5 * Y(0) * X(1) * Z(3) * Z(4) + Z(0) * Z(1) * X(2) * Z(3) * Z(4), nqubits=5) + circuit = Circuit(5) + circuit.add(gates.RX(0, theta=-np.pi / 2)) # = +1 + circuit.add(gates.H(1)) # = +1; q3,q4 in |0> -> exact = +0.5 + result = expectation_from_samples(circuit, hamiltonian, n_shots=100_000, grouping=method) + assert result == pytest.approx(hamiltonian.expectation(circuit), abs=0.01) + + def test_h2_hf_energy(): """Test HF energy of H2 molecule""" h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))]) From b9d934a1fb2208a25cf1733560817d57ddc280ab Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Fri, 24 Jul 2026 16:21:50 +0800 Subject: [PATCH 05/36] Update test code --- tests/test_measurement_util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_measurement_util.py b/tests/test_measurement_util.py index 529cdcb8..64994194 100644 --- a/tests/test_measurement_util.py +++ b/tests/test_measurement_util.py @@ -169,7 +169,7 @@ def test_lagrangian_subspace(): def test_sort_tau_terms(): # Using the example given in the function docstring test_terms = (["X0", "X2"], ["Z1", "X3", "Z4", "X5"], ["Z0", "Z2"], ["Z1"], ["Z3", "Z5"], ["Z4"]) - test_symplectic_form = [_pauli_to_symplectic(term, nqubits=6) for term in test_terms] + test_symplectic_form = np.array([_pauli_to_symplectic(term, nqubits=6) for term in test_terms]) result_symplectic = _sort_tau_terms(test_symplectic_form) result_pauli = [_symplectic_to_pauli(term) for term in result_symplectic] assert result_pauli == [["X0", "X2"], ["Z1"], ["Z0", "Z2"], ["Z3", "Z5"], ["Z4"], ["Z1", "X3", "Z4", "X5"]] From 125c14175e446af297359f61dff8fc3f9df9a1b7 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Fri, 24 Jul 2026 16:25:04 +0800 Subject: [PATCH 06/36] Loosen test bounds slightly --- tests/test_expectation_samples.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_expectation_samples.py b/tests/test_expectation_samples.py index 88808e66..251b4ee1 100644 --- a/tests/test_expectation_samples.py +++ b/tests/test_expectation_samples.py @@ -121,14 +121,14 @@ def test_measurement_grouping_extra_tests(): circuit.add(gates.H(0)) # = +1 circuit.add(gates.RX(1, theta=-np.pi / 2)) # = +1 -> exact = +0.5 result = expectation_from_samples(circuit, hamiltonian, n_shots=100_000, grouping=method) - assert result == pytest.approx(hamiltonian.expectation(circuit), abs=0.01) + assert result == pytest.approx(hamiltonian.expectation(circuit), abs=0.03) hamiltonian = SymbolicHamiltonian(0.5 * Y(0) * X(1) * Z(3) * Z(4) + Z(0) * Z(1) * X(2) * Z(3) * Z(4), nqubits=5) circuit = Circuit(5) circuit.add(gates.RX(0, theta=-np.pi / 2)) # = +1 circuit.add(gates.H(1)) # = +1; q3,q4 in |0> -> exact = +0.5 result = expectation_from_samples(circuit, hamiltonian, n_shots=100_000, grouping=method) - assert result == pytest.approx(hamiltonian.expectation(circuit), abs=0.01) + assert result == pytest.approx(hamiltonian.expectation(circuit), abs=0.03) def test_h2_hf_energy(): From 73395aa395351c7f082df4f616637a985409608d Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Wed, 29 Jul 2026 17:38:31 +0800 Subject: [PATCH 07/36] Minor improvements --- src/qibochem/measurement/util.py | 41 ++++++++++---------------------- tests/test_measurement_util.py | 6 ++--- 2 files changed, 16 insertions(+), 31 deletions(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index db2e8a72..a5fa2a4b 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -195,37 +195,22 @@ def _lagrangian_subspace(vector_space: np.ndarray) -> np.ndarray: def _sort_tau_terms(v_basis: np.ndarray) -> np.ndarray: - """ - Sorts v_basis s.t. the i'th term of basis vector i is NOT I, e.g. - [['X0', 'X2'], ['Z1', 'X3', 'Z4', 'X5'], ['Z0', 'Z2'], ['Z1'], ['Z3', 'Z5'], ['Z4']] - will return - [['X0', 'X2'], ['Z1'], ['Z0', 'Z2'], ['Z3', 'Z5'], ['Z4'], ['Z1', 'X3', 'Z4', 'X5']] - """ - # Check to see if sorting needed + """Sorts the rows of v_basis s.t. the (i, i) and (i, i+dim) entries are not 0, i.e. i'th basis vector i is NOT I""" dim = v_basis.shape[0] - if all(v_basis[i, i] or v_basis[i, i + dim] for i in range(dim)): - return v_basis - - # TODO: Can refactor this? - # Convert the basis set to strings for easier sorting - pauli_terms = [_symplectic_to_pauli(vector) for vector in v_basis] - dim = len(pauli_terms) - sorted_terms = {} - remaining = list(pauli_terms) - - while remaining: - qubit_terms = { - qubit: [term for term in remaining if any(_get_qubit(_op) == qubit for _op in term)] - for qubit in range(dim) - if qubit not in sorted_terms + while True: + # Sorting done + if all(v_basis[i, i] or v_basis[i, i + dim] for i in range(dim)): + break + # Sort unmatched qubits + unmatched_qubits = [i for i in range(dim) if not (v_basis[i, i] or v_basis[i, i + dim])] + matches_for_unmatched_qubits = { + i: [qubit for qubit in range(dim) if v_basis[i, qubit] or v_basis[i, qubit + dim]] for i in unmatched_qubits } # Preference: Qubits with fewest candidates (tie-break: min(qubit index)) - qubit = min(qubit_terms, key=lambda x: (len(qubit_terms[x]), x)) - selected_term = min(qubit_terms[qubit], key=len) - sorted_terms[qubit] = selected_term - remaining.remove(selected_term) - # Convert the strings back to symplectic vectors and return the whole array - return np.array([_pauli_to_symplectic(sorted_terms[_i], dim) for _i in range(dim)]) + row_to_swap = min(matches_for_unmatched_qubits, key=lambda x: (len(matches_for_unmatched_qubits[x]), x)) + target = min(matches_for_unmatched_qubits[row_to_swap]) + v_basis[[row_to_swap, target]] = v_basis[[target, row_to_swap]] + return v_basis def _get_sigma_terms(tau_terms: np.ndarray) -> tuple[np.ndarray, np.ndarray]: diff --git a/tests/test_measurement_util.py b/tests/test_measurement_util.py index 64994194..dd802ec6 100644 --- a/tests/test_measurement_util.py +++ b/tests/test_measurement_util.py @@ -168,15 +168,15 @@ def test_lagrangian_subspace(): def test_sort_tau_terms(): # Using the example given in the function docstring + nqubits = 6 test_terms = (["X0", "X2"], ["Z1", "X3", "Z4", "X5"], ["Z0", "Z2"], ["Z1"], ["Z3", "Z5"], ["Z4"]) test_symplectic_form = np.array([_pauli_to_symplectic(term, nqubits=6) for term in test_terms]) result_symplectic = _sort_tau_terms(test_symplectic_form) - result_pauli = [_symplectic_to_pauli(term) for term in result_symplectic] - assert result_pauli == [["X0", "X2"], ["Z1"], ["Z0", "Z2"], ["Z3", "Z5"], ["Z4"], ["Z1", "X3", "Z4", "X5"]] + assert all(result_symplectic[i, i] or result_symplectic[i, i + nqubits] for i in range(nqubits)) def test_get_sigma_terms(): - test_terms = [["X0", "X2"], ["Z1"], ["Z0", "Z2"], ["Z3", "Z5"], ["Z4"], ["Z1", "X3", "Z4", "X5"]] + test_terms = [["X0", "X2"], ["Z1"], ["Z0", "Z2"], ["Z1", "X3", "Z4", "X5"], ["Z4"], ["Z3", "Z5"]] test_symplectic_form = [_pauli_to_symplectic(term, nqubits=6) for term in test_terms] new_tau_terms, sigma_terms = _get_sigma_terms(test_symplectic_form) # Check new tau terms are still mutually orthogonal From d69726c337c19873d1afa8162296f14801bad91a Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Fri, 31 Jul 2026 16:54:24 +0800 Subject: [PATCH 08/36] Tentative fix to Gaussian elimination code --- src/qibochem/measurement/optimization.py | 1 + src/qibochem/measurement/util.py | 54 ++++++++++++++---------- 2 files changed, 32 insertions(+), 23 deletions(-) diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index 2d3ade2b..73e21f92 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -146,6 +146,7 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl qubit_map = dict(zip(sorted({q for pauli_op in x_result for q in pauli_op}), term_qubits)) phase_factors = [_phase_factor(v_basis[pauli_op]) for pauli_op in x_result] u_gates = _synthesise_circuit(v_basis) + u_gates += [gates.SWAP(i, j) for i, j in qubit_map.items() if i != j] mapping = { term: phase * prod(Z(qubit_map[_i]) for _i in soln) for term, phase, soln in zip(term_list, phase_factors, x_result) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index a5fa2a4b..6703e1af 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -325,34 +325,39 @@ def _col_reduce_x_matrix(stabiliser_matrix: np.ndarray) -> list[gates.Gate]: list[gates.Gate]: List of CNOT/SWAP gates to be added to the circuit """ gates_list = [] - dim, _dim_space = stabiliser_matrix.shape - dim_space = _dim_space // 2 + dim, dim_space = stabiliser_matrix.shape + dim_space = dim_space // 2 + pivot_col = 0 # Paper used row reduction, but should be column reduction in our context - for _i in range(dim_space): - if _i >= dim: + for row in range(dim): + if pivot_col >= dim_space: break - # Get columns with row _i != 0 - nonzero_cols = np.nonzero(stabiliser_matrix[_i, :dim_space])[0] - - # Always take the first nonzero row to sort - _col = [_j for _j in nonzero_cols if _j >= _i][0] - if _i not in nonzero_cols: - stabiliser_matrix[:, [_i, _col, _i + dim_space, _col + dim_space]] = stabiliser_matrix[ - :, [_col, _i, _col + dim_space, _i + dim_space] + # Get columns at row i with 1 + nonzero_cols = np.where(stabiliser_matrix[row, pivot_col:dim_space] == 1)[0] + + col = pivot_col + nonzero_cols[0] + + # Move pivot column of X matrix into position + if col != pivot_col: + stabiliser_matrix[:, [pivot_col, col, pivot_col + dim_space, col + dim_space]] = stabiliser_matrix[ + :, [col, pivot_col, col + dim_space, pivot_col + dim_space] ] - gates_list.append(gates.SWAP(_i, _col)) - nonzero_cols = np.nonzero(stabiliser_matrix[_i, :dim_space])[0] + gates_list.append(gates.SWAP(col, pivot_col)) + + # Eliminate other 1's in the present row + nonzero_cols = np.where(stabiliser_matrix[row, :dim_space] == 1)[0] + nonzero_cols = nonzero_cols[nonzero_cols != pivot_col] + # Remove all nonzero entries on row _i using CNOT gates - for _col in nonzero_cols: # Ignore first entry of nonzero_cols since effectively should be 0 now - if _col != _i: - # Add j^th column to i^th column - stabiliser_matrix[:, _col] += stabiliser_matrix[:, _i] - # Add (i+dim_space)^th column to (j+dim_space)^th column - # RHS of stabiliser matrix should be 0 matrix, so I think can ignore...? - stabiliser_matrix[:, dim_space + _i] += stabiliser_matrix[:, _col + dim_space] - stabiliser_matrix %= 2 - gates_list.append(gates.CNOT(_i, _col)) + for col in nonzero_cols: + # X matrix: Add pivot column to column with 1 + stabiliser_matrix[:, col] += stabiliser_matrix[:, pivot_col] + # Z matrix: Add (column with 1)^th column to pivot column + stabiliser_matrix[:, pivot_col + dim_space] += stabiliser_matrix[:, col + dim_space] + stabiliser_matrix %= 2 + gates_list.append(gates.CNOT(col, pivot_col)) + pivot_col += 1 return gates_list @@ -397,10 +402,13 @@ def _synthesise_circuit(v_basis: np.ndarray) -> list[gates.Gate]: rotation_gates = [] # 1. Apply H gates to transform 'X matrix' to full rank rotation_gates += _make_x_matrix_full_rank(stabiliser_matrix) + print("Matrix:\n", stabiliser_matrix) # 2. Row-reduce 'X matrix' to I using CNOT/SWAP gates rotation_gates += _col_reduce_x_matrix(stabiliser_matrix) + print("Matrix:\n", stabiliser_matrix) # 3. Remove all non-zero entries on 'Z matrix' using S and CZ gates rotation_gates += _zero_z_matrix(stabiliser_matrix) + print("Matrix:\n", stabiliser_matrix) # 4. Apply H to each qubit to swap the 'X' and 'Z' matrices rotation_gates += [gates.H(_i) for _i in range(n_qubits)] return rotation_gates From 2c7f513e8f24c00d15c019d58bfb9526660d9f63 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 3 Aug 2026 11:53:53 +0800 Subject: [PATCH 09/36] Change symplectic arrays dtype: int -> np.uint8 --- src/qibochem/measurement/util.py | 49 ++++++++++++++++++----- tests/test_measurement_util.py | 69 +++++++++++++++++++++----------- 2 files changed, 85 insertions(+), 33 deletions(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index 6703e1af..42db56c6 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -87,7 +87,9 @@ def _pauli_to_symplectic(pauli_string: list[str], nqubits: int) -> np.ndarray: pauli_ops = {_get_qubit(pauli_op): pauli_op[0] for pauli_op in pauli_string} # Pauli operator for each qubit # Convert to the symplectic vector sym_vector = np.reshape( - np.array([PAULI_BINARY[pauli_ops.get(_i, "I")] for _i in range(nqubits)]), shape=2 * nqubits, order="F" + np.array([PAULI_BINARY[pauli_ops.get(i, "I")] for i in range(nqubits)], dtype=np.uint8), + shape=2 * nqubits, + order="F", ) return sym_vector @@ -115,10 +117,36 @@ def _symplectic_inner_product(u: np.ndarray, v: np.ndarray) -> int: def _binary_gaussian_elimination(vector_space: np.ndarray) -> np.ndarray: """ - Carries out Gaussian elimination on a binary vector_space to obtain a basis for vector_space. Reduces and returns - vector_space to its (unique) reduced row echelon form, and removes any zero rows as well + Carries out Gaussian elimination on a binary vector_space to obtain a basis for vector_space. Reduces vector_space + in-place to its (unique) reduced row echelon form, and removes any zero rows as well """ - cp_vector_space = np.array(vector_space) + cp_vector_space = np.array(vector_space, dtype=np.uint8) + # rows, cols = vector_space.shape + + # pivot_row = 0 + # for col in range(cols): + # # Find a pivot row with a 1 in current column. + # pivot_candidates = np.where(vector_space[pivot_row:, col] == 1)[0] + # if pivot_candidates.size == 0: + # continue + + # row = pivot_row + pivot_candidates[0] + + # # Swap current row with pivot row if needed. + # if pivot_row != row: + # vector_space[[row, pivot_row]] = vector_space[[pivot_row, row]] + + # # Eliminate all other rows + # rows_to_reduce = np.where(vector_space[:, col] == 1)[0] + # rows_to_reduce = rows_to_reduce[rows_to_reduce != row] + + # # In GF(2), elimination is XOR with the pivot row. + # vector_space[rows_to_reduce] += vector_space[row] + # vector_space %= 2 + + # row += 1 + # if row == rows: + # break dim = vector_space.shape[0] # Swap the rows in the vector space to get its row echelon form @@ -144,8 +172,11 @@ def _binary_gaussian_elimination(vector_space: np.ndarray) -> np.ndarray: # Remove all zero rows from the obtained basis zero_vector_indices = np.all(cp_vector_space == 0, axis=1) cp_vector_space = cp_vector_space[~zero_vector_indices] + return cp_vector_space + # return vector_space + def _binary_nullspace(binary_matrix: np.ndarray) -> np.ndarray: """Finds the nullspace of a binary_matrix, i.e. x s.t. Ax = 0""" @@ -222,7 +253,7 @@ def _get_sigma_terms(tau_terms: np.ndarray) -> tuple[np.ndarray, np.ndarray]: sigma_terms = [] dim = tau_terms[0].shape[0] // 2 # Make a copy of the original basis set for orthogonalization - new_tau_terms = np.array(tau_terms) + new_tau_terms = np.array(tau_terms, dtype=np.uint8) # Iterate over the original tau_i to make changes to new_tau_i for _i in range(dim): tau_i = new_tau_terms[_i] @@ -239,8 +270,9 @@ def _get_sigma_terms(tau_terms: np.ndarray) -> tuple[np.ndarray, np.ndarray]: _symplectic_inner_product(new_tau_terms[_j], sigma_i) * tau_i if _j != _i else np.zeros(2 * dim) # symplectic_inner_product(new_tau_terms[_j], sigma_i) * tau_i if _j > _i else np.zeros(2 * dim) for _j in range(dim) - ] - ).astype(int) + ], + dtype=np.uint8, + ) new_tau_terms = new_tau_terms % 2 return new_tau_terms, np.array(sigma_terms) @@ -402,13 +434,10 @@ def _synthesise_circuit(v_basis: np.ndarray) -> list[gates.Gate]: rotation_gates = [] # 1. Apply H gates to transform 'X matrix' to full rank rotation_gates += _make_x_matrix_full_rank(stabiliser_matrix) - print("Matrix:\n", stabiliser_matrix) # 2. Row-reduce 'X matrix' to I using CNOT/SWAP gates rotation_gates += _col_reduce_x_matrix(stabiliser_matrix) - print("Matrix:\n", stabiliser_matrix) # 3. Remove all non-zero entries on 'Z matrix' using S and CZ gates rotation_gates += _zero_z_matrix(stabiliser_matrix) - print("Matrix:\n", stabiliser_matrix) # 4. Apply H to each qubit to swap the 'X' and 'Z' matrices rotation_gates += [gates.H(_i) for _i in range(n_qubits)] return rotation_gates diff --git a/tests/test_measurement_util.py b/tests/test_measurement_util.py index dd802ec6..936d3337 100644 --- a/tests/test_measurement_util.py +++ b/tests/test_measurement_util.py @@ -57,9 +57,9 @@ def test_group_commuting_terms(term_list, qwc_expected, gc_expected): # "pauli_string,n_qubits,expected", "function_args,expected", [ - ({"pauli_string": ["X0", "Y1", "Z2"], "nqubits": 4}, np.array([1, 1, 0, 0, 0, 1, 1, 0])), - ({"pauli_string": ["Z1", "X3"], "nqubits": 4}, np.array([0, 0, 0, 1, 0, 1, 0, 0])), - ({"pauli_string": [], "nqubits": 4}, np.array([0, 0, 0, 0, 0, 0, 0, 0])), + ({"pauli_string": ["X0", "Y1", "Z2"], "nqubits": 4}, np.array([1, 1, 0, 0, 0, 1, 1, 0], dtype=np.uint8)), + ({"pauli_string": ["Z1", "X3"], "nqubits": 4}, np.array([0, 0, 0, 1, 0, 1, 0, 0], dtype=np.uint8)), + ({"pauli_string": [], "nqubits": 4}, np.array([0, 0, 0, 0, 0, 0, 0, 0], dtype=np.uint8)), ], ) def test_pauli_to_symplectic(function_args, expected): @@ -70,8 +70,8 @@ def test_pauli_to_symplectic(function_args, expected): @pytest.mark.parametrize( "function_args,expected", [ - ({"symplectic_vector": np.array([1, 1, 0, 0, 0, 1, 1, 0])}, ["X0", "Y1", "Z2"]), - ({"symplectic_vector": np.array([0, 1, 0, 1, 0, 1, 0, 0])}, ["Y1", "X3"]), + ({"symplectic_vector": np.array([1, 1, 0, 0, 0, 1, 1, 0], dtype=np.uint8)}, ["X0", "Y1", "Z2"]), + ({"symplectic_vector": np.array([0, 1, 0, 1, 0, 1, 0, 0], dtype=np.uint8)}, ["Y1", "X3"]), ], ) def test_symplectic_to_pauli(function_args, expected): @@ -82,8 +82,8 @@ def test_symplectic_to_pauli(function_args, expected): @pytest.mark.parametrize( "u,v", [ - (np.array([1, 1, 0, 0, 0, 1, 1, 0]), np.array([1, 1, 0, 0, 0, 1, 1, 0])), - (np.array([1, 0, 0, 0, 1, 1, 1, 1]), np.array([1, 1, 0, 0, 0, 1, 1, 0])), + (np.array([1, 1, 0, 0, 0, 1, 1, 0], dtype=np.uint8), np.array([1, 1, 0, 0, 0, 1, 1, 0], dtype=np.uint8)), + (np.array([1, 0, 0, 0, 1, 1, 1, 1], dtype=np.uint8), np.array([1, 1, 0, 0, 0, 1, 1, 0], dtype=np.uint8)), ], ) def test_symplectic_inner_product(u, v): @@ -91,12 +91,12 @@ def test_symplectic_inner_product(u, v): dim = u.shape[0] // 2 j_matrix = np.concatenate( ( - np.concatenate((np.zeros((dim, dim)), np.identity(dim, dtype=int)), axis=1), - np.concatenate((np.identity(dim, dtype=int), np.zeros((dim, dim))), axis=1), + np.concatenate((np.zeros((dim, dim), dtype=np.uint8), np.identity(dim, dtype=np.uint8)), axis=1), + np.concatenate((np.identity(dim, dtype=np.uint8), np.zeros((dim, dim), dtype=np.uint8)), axis=1), ), axis=0, ) - assert _symplectic_inner_product(u, v) == (np.dot(u, np.dot(j_matrix, v)).astype(int) % 2) + assert _symplectic_inner_product(u, v) == (np.dot(u, np.dot(j_matrix, v)) % 2) @pytest.mark.parametrize( @@ -104,30 +104,49 @@ def test_symplectic_inner_product(u, v): [ ( np.array( - [[0, 1, 1, 0, 0, 0], [1, 1, 0, 0, 1, 1], [0, 1, 1, 0, 0, 0], [1, 1, 0, 0, 1, 1], [0, 0, 1, 0, 1, 1]] + [[0, 1, 1, 0, 0, 0], [1, 1, 0, 0, 1, 1], [0, 1, 1, 0, 0, 0], [1, 1, 0, 0, 1, 1], [0, 0, 1, 0, 1, 1]], + dtype=np.uint8, ), - np.array([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 1, 1], [0, 0, 1, 0, 1, 1]]), + np.array([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 1, 1], [0, 0, 1, 0, 1, 1]], dtype=np.uint8), ), ( np.array( - [[1, 1, 1, 1, 0, 1, 1, 0], [1, 1, 1, 1, 1, 0, 0, 1], [1, 1, 1, 1, 0, 0, 1, 1], [1, 1, 1, 1, 1, 1, 0, 0]] + [ + [1, 1, 1, 1, 0, 1, 1, 0], + [1, 1, 1, 1, 1, 0, 0, 1], + [1, 1, 1, 1, 0, 0, 1, 1], + [1, 1, 1, 1, 1, 1, 0, 0], + ], + dtype=np.uint8, ), - np.array([[1, 1, 1, 1, 0, 0, 1, 1], [0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 0, 1, 0, 1]]), + np.array([[1, 1, 1, 1, 0, 0, 1, 1], [0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 0, 1, 0, 1]], dtype=np.uint8), ), ( - np.array([[0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 1, 1, 1, 1], [0, 0, 0, 0, 0, 1, 0, 1]]), - np.array([[0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 0, 1, 0, 1]]), + np.array([[0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 1, 1, 1, 1], [0, 0, 0, 0, 0, 1, 0, 1]], dtype=np.uint8), + np.array([[0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 0, 1, 0, 1]], dtype=np.uint8), ), ( - np.array([[0, 1, 0, 1, 0, 1], [0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 1, 0]]), - np.array([[0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 1, 0]]), + np.array([[0, 1, 0, 1, 0, 1], [0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 1, 0]], dtype=np.uint8), + np.array([[0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 1, 0]], dtype=np.uint8), ), ( np.array( - [[1, 1, 1, 1, 0, 0, 0, 0], [1, 1, 1, 1, 0, 0, 1, 1], [1, 1, 1, 1, 0, 1, 0, 1], [1, 1, 1, 1, 1, 0, 0, 1]] + [ + [1, 1, 1, 1, 0, 0, 0, 0], + [1, 1, 1, 1, 0, 0, 1, 1], + [1, 1, 1, 1, 0, 1, 0, 1], + [1, 1, 1, 1, 1, 0, 0, 1], + ], + dtype=np.uint8, ), np.array( - [[1, 1, 1, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 1], [0, 0, 0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 0, 0, 1, 1]] + [ + [1, 1, 1, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 1], + [0, 0, 0, 0, 0, 1, 0, 1], + [0, 0, 0, 0, 0, 0, 1, 1], + ], + dtype=np.uint8, ), ), ], @@ -139,9 +158,13 @@ def test_binary_gaussian_elimination(test, result): def test_binary_nullspace(): - test_space = np.array([[1, 1, 1, 1, 0, 0, 1, 1], [0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 0, 1, 0, 1]]) + test_space = np.array( + [[1, 1, 1, 1, 0, 0, 1, 1], [0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 0, 1, 0, 1]], dtype=np.uint8 + ) nullspace = _binary_nullspace(test_space) - assert all(np.allclose((test_space @ vector) % 2, np.zeros(test_space.shape[0])) for vector in nullspace) + assert all( + np.allclose((test_space @ vector) % 2, np.zeros(test_space.shape[0], dtype=np.uint8)) for vector in nullspace + ) def test_lagrangian_subspace(): @@ -154,7 +177,7 @@ def test_lagrangian_subspace(): [0, 0, 0, 1, 0, 1, 0, 1], [0, 0, 0, 0, 1, 1, 1, 1], ], - dtype=int, + dtype=np.uint8, ) subspace = _lagrangian_subspace(test_space) # Vectors in subspace should all be symplectically orthogonal to each other From 7b98afaa613f97c406d0850dc4908b59cb913fdf Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 3 Aug 2026 15:41:12 +0800 Subject: [PATCH 10/36] Refactor some symplectic array-related functions --- src/qibochem/measurement/optimization.py | 2 +- src/qibochem/measurement/util.py | 111 +++++++++-------------- 2 files changed, 43 insertions(+), 70 deletions(-) diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index 73e21f92..c0ce366e 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -126,7 +126,7 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl # Otherwise, expression is a sum of terms term_list = [_term_to_string(term) for term in expression.args if _term_to_string(term)[0] in ("X", "Y", "Z")] term_qubits = sorted({_get_qubit(op) for term in term_list for op in term.split()}) - v_subspace = np.array([_pauli_to_symplectic(terms.split(), nqubits) for terms in term_list]) + v_subspace = np.array([_pauli_to_symplectic(terms.split(), nqubits) for terms in term_list], dtype=np.uint8) v_basis = _binary_gaussian_elimination(v_subspace) dim_v = v_basis.shape[0] diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index 42db56c6..c81c93c5 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -5,6 +5,7 @@ import networkx as nx import numpy as np from qibo import gates +from qibo.config import raise_error # Mapping of Pauli operators to a symplectic (binary) representation, folowing the convention of (X|Z) PAULI_BINARY = {"I": (0, 0), "X": (1, 0), "Y": (1, 1), "Z": (0, 1)} @@ -117,65 +118,41 @@ def _symplectic_inner_product(u: np.ndarray, v: np.ndarray) -> int: def _binary_gaussian_elimination(vector_space: np.ndarray) -> np.ndarray: """ - Carries out Gaussian elimination on a binary vector_space to obtain a basis for vector_space. Reduces vector_space - in-place to its (unique) reduced row echelon form, and removes any zero rows as well + Performs Gaussian elimination on a binary vector_space. Returns the (unique) reduced row echelon form, and removes + any zero rows as well """ - cp_vector_space = np.array(vector_space, dtype=np.uint8) - # rows, cols = vector_space.shape - - # pivot_row = 0 - # for col in range(cols): - # # Find a pivot row with a 1 in current column. - # pivot_candidates = np.where(vector_space[pivot_row:, col] == 1)[0] - # if pivot_candidates.size == 0: - # continue - - # row = pivot_row + pivot_candidates[0] + vector_space = np.array(vector_space) # Create a copy for returning + rows, cols = vector_space.shape + + pivot_row = 0 + for col in range(cols): + # Find a pivot row with a 1 in current column. + pivot_candidates = np.where(vector_space[pivot_row:, col] == 1)[0] + if pivot_candidates.size == 0: + continue - # # Swap current row with pivot row if needed. - # if pivot_row != row: - # vector_space[[row, pivot_row]] = vector_space[[pivot_row, row]] + row = pivot_row + pivot_candidates[0] - # # Eliminate all other rows - # rows_to_reduce = np.where(vector_space[:, col] == 1)[0] - # rows_to_reduce = rows_to_reduce[rows_to_reduce != row] + # Swap current row with pivot row if needed. + if pivot_row != row: + vector_space[[row, pivot_row]] = vector_space[[pivot_row, row]] - # # In GF(2), elimination is XOR with the pivot row. - # vector_space[rows_to_reduce] += vector_space[row] - # vector_space %= 2 + # Eliminate all other rows + rows_to_reduce = np.where(vector_space[:, col] == 1)[0] + rows_to_reduce = rows_to_reduce[rows_to_reduce != pivot_row] - # row += 1 - # if row == rows: - # break + # In GF(2), elimination is XOR with the pivot row. + vector_space[rows_to_reduce] += vector_space[pivot_row] + vector_space %= 2 - dim = vector_space.shape[0] - # Swap the rows in the vector space to get its row echelon form - for _i in range(dim): - subspace_to_sort = cp_vector_space[_i:, :] - if not np.any(subspace_to_sort): + pivot_row += 1 + if pivot_row == rows: break - # Always take the first nonzero column to sort - nonzero_cols = np.nonzero(np.any(subspace_to_sort, axis=0))[0] - _col = nonzero_cols[0] - - col_indices = subspace_to_sort[:, _col].argsort()[::-1] - subspace_to_sort[:, :] = subspace_to_sort[col_indices] - - # Other than row _i, find which rows have 1 in column _i - rows_to_reduce = [_j for _j, _row in enumerate(cp_vector_space) if _j != _i and cp_vector_space[_j, _col] == 1] - - # Add row _i to each of the rows with 1 in the same column - cp_vector_space[[rows_to_reduce]] += cp_vector_space[_i] - cp_vector_space %= 2 - # Remove all zero rows from the obtained basis - zero_vector_indices = np.all(cp_vector_space == 0, axis=1) - cp_vector_space = cp_vector_space[~zero_vector_indices] - - return cp_vector_space - - # return vector_space + zero_vector_indices = np.all(vector_space == 0, axis=1) + vector_space = vector_space[~zero_vector_indices] + return vector_space def _binary_nullspace(binary_matrix: np.ndarray) -> np.ndarray: @@ -228,10 +205,7 @@ def _lagrangian_subspace(vector_space: np.ndarray) -> np.ndarray: def _sort_tau_terms(v_basis: np.ndarray) -> np.ndarray: """Sorts the rows of v_basis s.t. the (i, i) and (i, i+dim) entries are not 0, i.e. i'th basis vector i is NOT I""" dim = v_basis.shape[0] - while True: - # Sorting done - if all(v_basis[i, i] or v_basis[i, i + dim] for i in range(dim)): - break + while not all(v_basis[i, i] or v_basis[i, i + dim] for i in range(dim)): # Sort unmatched qubits unmatched_qubits = [i for i in range(dim) if not (v_basis[i, i] or v_basis[i, i + dim])] matches_for_unmatched_qubits = { @@ -275,7 +249,7 @@ def _get_sigma_terms(tau_terms: np.ndarray) -> tuple[np.ndarray, np.ndarray]: ) new_tau_terms = new_tau_terms % 2 - return new_tau_terms, np.array(sigma_terms) + return new_tau_terms, np.array(sigma_terms, dtype=np.uint8) def _solve_linear_system(binary_matrix: np.ndarray, vector: np.ndarray) -> list[np.ndarray]: @@ -284,7 +258,7 @@ def _solve_linear_system(binary_matrix: np.ndarray, vector: np.ndarray) -> list[ aug_matrix = np.concatenate((binary_matrix, vector), axis=0).T rref_aug_matrix = _binary_gaussian_elimination(aug_matrix) # Get non-zero entries in each column on RHS of rref_aug_matrix => Solution for respective vector in b - return [np.nonzero(rref_aug_matrix[:, binary_matrix.shape[0] + _i])[0].tolist() for _i in range(vector.shape[0])] + return [np.nonzero(rref_aug_matrix[:, binary_matrix.shape[0] + i])[0].tolist() for i in range(vector.shape[0])] def _single_qubit_phase_factor(pauli_ops: list[np.ndarray]) -> complex: @@ -332,20 +306,19 @@ def _make_x_matrix_full_rank(stabiliser_matrix: np.ndarray) -> list[gates.Gate]: dim_space = stabiliser_matrix.shape[1] // 2 x_matrix = stabiliser_matrix[:, :dim_space] z_matrix = stabiliser_matrix[:, dim_space:] + # Need to find full rank submatrix in Z matrix for each of the zero rows in the X matrix - zero_row_indices = [_i for _i, is_zero in enumerate(np.all(x_matrix == 0, axis=1)) if is_zero] - # Only need to do anything if there are zero rows in the X matrix - while zero_row_indices: - nonzero_cols_by_row = {row: list(np.nonzero(z_matrix[row, :])[0]) for row in zero_row_indices} - # See if there are any single-element lists in the values of nonzero_cols_by_row - no_choice_rows = [row for row, possible_cols in nonzero_cols_by_row.items() if len(possible_cols) == 1] - chosen_row = no_choice_rows[0] if no_choice_rows else zero_row_indices[0] - chosen_qubit = nonzero_cols_by_row[chosen_row][0] - stabiliser_matrix[:, [chosen_qubit, chosen_qubit + dim_space]] = stabiliser_matrix[ - :, [chosen_qubit + dim_space, chosen_qubit] - ] - gates_list.append(gates.H(chosen_qubit)) - zero_row_indices.remove(chosen_row) + zero_row_indices = np.where(np.all(x_matrix == 0, axis=1))[0] + prev_swap = [] # To avoid swapping the same column twice + while zero_row_indices.size > 0: + # Select the first possible column for the first zero row + for qubit in np.nonzero(z_matrix[zero_row_indices[0], :])[0]: + if qubit not in prev_swap: + stabiliser_matrix[:, [qubit, qubit + dim_space]] = stabiliser_matrix[:, [qubit + dim_space, qubit]] + gates_list.append(gates.H(qubit)) + prev_swap.append(qubit) + break + zero_row_indices = np.where(np.all(x_matrix == 0, axis=1))[0] return gates_list From ec2eee1b1087a4080ce3b6c986a39d6720690d19 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 3 Aug 2026 16:05:59 +0800 Subject: [PATCH 11/36] Cleanup tests (in-progress) --- tests/test_expectation_samples.py | 112 +++++++++++++++--------------- tests/test_measurement_util.py | 6 +- 2 files changed, 61 insertions(+), 57 deletions(-) diff --git a/tests/test_expectation_samples.py b/tests/test_expectation_samples.py index 251b4ee1..af209901 100644 --- a/tests/test_expectation_samples.py +++ b/tests/test_expectation_samples.py @@ -32,20 +32,20 @@ def test_pauli_term_measurement_expectation(term, frequencies, qubit_map, expect assert result == expected -@pytest.mark.parametrize( - "terms,gates_to_add", - [ - (Z(0), [gates.X(0)]), - (Z(0) * Z(1), [gates.X(0)]), - (X(0), [gates.H(0)]), - ], -) -def test_expectation_from_samples(terms, gates_to_add): - hamiltonian = SymbolicHamiltonian(terms, nqubits=2) - circuit = Circuit(2) - circuit.add(gates_to_add) - result = expectation_from_samples(circuit, hamiltonian) - assert result == pytest.approx(expected := hamiltonian.expectation(circuit)), f"{result} != {expected}" +# @pytest.mark.parametrize( +# "terms,gates_to_add", +# [ +# (Z(0), [gates.X(0)]), +# (Z(0) * Z(1), [gates.X(0)]), +# (X(0), [gates.H(0)]), +# ], +# ) +# def test_expectation_from_samples(terms, gates_to_add): +# hamiltonian = SymbolicHamiltonian(terms, nqubits=2) +# circuit = Circuit(2) +# circuit.add(gates_to_add) +# result = expectation_from_samples(circuit, hamiltonian) +# assert result == pytest.approx(expected := hamiltonian.expectation(circuit)), f"{result} != {expected}" def test_measurement_basis_rotations_error(): @@ -82,53 +82,55 @@ def test_expectation_invalid_shot_allocation(): ) +@pytest.mark.parametrize("grouping", ["qwc", "gc", "gc2"]) @pytest.mark.parametrize( - "hamiltonian", + "terms", [ - SymbolicHamiltonian(Z(2)), - SymbolicHamiltonian(0.2 * X(0) + Y(2) + 13.0), - SymbolicHamiltonian(Z(0) + X(0) * Y(1) + Z(0) * Y(2)), - SymbolicHamiltonian(Y(0) + Z(1) + X(0) * Z(2)), - SymbolicHamiltonian( - 0.1 * X(0) * X(1) * Y(2) + 0.2 * X(0) * Y(1) * Y(2) + 0.3 * Y(0) * X(1) * X(2) - 3.14 * Y(0) * Y(1) * X(2) - ), + Z(2), + 0.2 * X(0) + Y(2) + 13.0, + Z(0) + X(0) * Y(1) + Z(0) * Y(2), + Y(0) + Z(1) + X(0) * Z(2), + 0.1 * X(0) * X(1) * Y(2) + 0.2 * X(0) * Y(1) * Y(2) + 0.3 * Y(0) * X(1) * X(2) - 3.14 * Y(0) * Y(1) * X(2), ], ) -def test_measurement_grouping_functionality(hamiltonian): +def test_measurement_grouping_functionality(grouping, terms): """Small scale tests of commuting measurements functionality""" - n_qubits = 3 - circuit = Circuit(n_qubits) - circuit.add(gates.RX(_i, 0.1 * _i) for _i in range(n_qubits)) - circuit.add(gates.CNOT(_i, _i + 1) for _i in range(n_qubits - 1)) - circuit.add(gates.RZ(_i, 0.2 * _i) for _i in range(n_qubits)) + nqubits = 3 + circuit = Circuit(nqubits) + circuit.add(gates.RX(_i, 0.1 * _i) for _i in range(nqubits)) + circuit.add(gates.CNOT(_i, _i + 1) for _i in range(nqubits - 1)) + circuit.add(gates.RZ(_i, 0.2 * _i) for _i in range(nqubits)) + hamiltonian = SymbolicHamiltonian(terms, nqubits=nqubits) expected = hamiltonian.expectation(circuit) - n_shots = 10000 - for grouping in ("qwc", "gc", "gc2"): - test = expectation_from_samples( - circuit, - hamiltonian, - n_shots=n_shots, - grouping=grouping, - ) - assert test == pytest.approx(expected, abs=0.08) - - -def test_measurement_grouping_extra_tests(): - """Additional tests for some previously undiscovered bugs""" - for method in ("gc", "gc2"): - hamiltonian = SymbolicHamiltonian(0.5 * X(0) * Y(1) + Z(0) * Z(1) * Z(2), nqubits=3) - circuit = Circuit(3) - circuit.add(gates.H(0)) # = +1 - circuit.add(gates.RX(1, theta=-np.pi / 2)) # = +1 -> exact = +0.5 - result = expectation_from_samples(circuit, hamiltonian, n_shots=100_000, grouping=method) - assert result == pytest.approx(hamiltonian.expectation(circuit), abs=0.03) - - hamiltonian = SymbolicHamiltonian(0.5 * Y(0) * X(1) * Z(3) * Z(4) + Z(0) * Z(1) * X(2) * Z(3) * Z(4), nqubits=5) - circuit = Circuit(5) - circuit.add(gates.RX(0, theta=-np.pi / 2)) # = +1 - circuit.add(gates.H(1)) # = +1; q3,q4 in |0> -> exact = +0.5 - result = expectation_from_samples(circuit, hamiltonian, n_shots=100_000, grouping=method) - assert result == pytest.approx(hamiltonian.expectation(circuit), abs=0.03) + test = expectation_from_samples( + circuit, + hamiltonian, + n_shots=100_000, + grouping=grouping, + ) + assert test == pytest.approx(expected, abs=0.05) + + +@pytest.mark.parametrize("grouping", ["gc", "gc2"]) +@pytest.mark.parametrize( + "terms,nqubits,gates_to_add", + [ + (0.5 * X(0) * Y(1) + Z(0) * Z(1) * Z(2), 3, (gates.H(0), gates.RX(1, theta=-np.pi / 2))), + ( + 0.5 * Y(0) * X(1) * Z(3) * Z(4) + Z(0) * Z(1) * X(2) * Z(3) * Z(4), + 5, + (gates.RX(0, theta=-np.pi / 2), gates.H(1)), + ), + (Y(1) * Y(2) + X(0) * X(1) * Z(2), 3, (gates.H(0), gates.H(1))), + ], +) +def test_measurement_grouping_extra_tests(grouping, terms, nqubits, gates_to_add): + """Additional tests for generally commuting terms""" + hamiltonian = SymbolicHamiltonian(terms, nqubits=nqubits) + circuit = Circuit(nqubits) + circuit.add(gates_to_add) + result = expectation_from_samples(circuit, hamiltonian, n_shots=100_000, grouping=grouping) + assert result == pytest.approx(hamiltonian.expectation(circuit), abs=0.03) def test_h2_hf_energy(): diff --git a/tests/test_measurement_util.py b/tests/test_measurement_util.py index 936d3337..5763c31a 100644 --- a/tests/test_measurement_util.py +++ b/tests/test_measurement_util.py @@ -254,7 +254,8 @@ def test_col_reduce_x_matrix(): [ [1, 1, 0, 0, 0, 1, 0, 1], [0, 1, 0, 0, 1, 0, 1, 0], - ] + ], + dtype=np.uint8, ) gates_list = _col_reduce_x_matrix(stabiliser_matrix) # Single column operation, should have only CNOT gate @@ -268,7 +269,8 @@ def test_zero_z_matrix(): [0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0], - ] + ], + dtype=np.uint8, ) gates_list = _zero_z_matrix(stabiliser_matrix) # Single column operation, should have only CNOT gate From a5f7260ea82cca248be843e0755323c65477bd56 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 3 Aug 2026 16:15:47 +0800 Subject: [PATCH 12/36] Finish cleaning up tests for sample stuff --- tests/test_expectation_samples.py | 56 ++++++++++++------------------- 1 file changed, 21 insertions(+), 35 deletions(-) diff --git a/tests/test_expectation_samples.py b/tests/test_expectation_samples.py index af209901..effe6495 100644 --- a/tests/test_expectation_samples.py +++ b/tests/test_expectation_samples.py @@ -32,22 +32,6 @@ def test_pauli_term_measurement_expectation(term, frequencies, qubit_map, expect assert result == expected -# @pytest.mark.parametrize( -# "terms,gates_to_add", -# [ -# (Z(0), [gates.X(0)]), -# (Z(0) * Z(1), [gates.X(0)]), -# (X(0), [gates.H(0)]), -# ], -# ) -# def test_expectation_from_samples(terms, gates_to_add): -# hamiltonian = SymbolicHamiltonian(terms, nqubits=2) -# circuit = Circuit(2) -# circuit.add(gates_to_add) -# result = expectation_from_samples(circuit, hamiltonian) -# assert result == pytest.approx(expected := hamiltonian.expectation(circuit)), f"{result} != {expected}" - - def test_measurement_basis_rotations_error(): """If unknown measurement grouping scheme used""" hamiltonian = SymbolicHamiltonian(Z(0) + X(0)) @@ -158,18 +142,19 @@ def test_h2_hf_energy(): @pytest.mark.parametrize( - "hamiltonian,grouping,expected_means,expected_variances", + "terms,grouping,expected_means,expected_variances", [ - (SymbolicHamiltonian(X(0), nqubits=2), None, [1.0], [0.0]), - (SymbolicHamiltonian(X(0) + Z(0), nqubits=2), None, [1.0, 0.0], [0.0, 0.0]), - (SymbolicHamiltonian(Z(0) + X(0) * Z(1), nqubits=2), "qwc", [-1.0, 0.0], [0.0, 0.0]), + (X(0), None, [1.0], [0.0]), + (X(0) + Z(0), None, [1.0, 0.0], [0.0, 0.0]), + (Z(0) + X(0) * Z(1), "qwc", [-1.0, 0.0], [0.0, 0.0]), ], ) -def test_sample_statistics(hamiltonian, grouping, expected_means, expected_variances): +def test_sample_statistics(terms, grouping, expected_means, expected_variances): circuit = Circuit(2) circuit.add(gates.H(0)) circuit.add(gates.X(1)) - n_trial_shots = 20000 + n_trial_shots = 20_000 + hamiltonian = SymbolicHamiltonian(terms, nqubits=2) grouped_terms = _measurement_basis_rotations(hamiltonian, grouping) sample_means, sample_variances = sample_statistics(circuit, grouped_terms, n_shots=n_trial_shots) assert sample_means == pytest.approx(expected_means, abs=0.08) @@ -177,24 +162,25 @@ def test_sample_statistics(hamiltonian, grouping, expected_means, expected_varia @pytest.mark.parametrize( - "hamiltonian,grouping", + "terms,grouping", [ - (SymbolicHamiltonian(0.2 * X(0) + Y(2) + 13.0), None), - (SymbolicHamiltonian(0.2 * X(0) + Y(2) + 13.0), "qwc"), - (SymbolicHamiltonian(Z(0) + X(0) * Y(1) + Z(0) * Y(2)), None), - (SymbolicHamiltonian(Y(0) + Z(1) + X(0) * Z(2)), "qwc"), - (SymbolicHamiltonian(Y(0) + Z(1) + X(0) * Z(2)), "gc"), + (0.2 * X(0) + Y(2) + 13.0, None), + (0.2 * X(0) + Y(2) + 13.0, "qwc"), + (Z(0) + X(0) * Y(1) + Z(0) * Y(2), None), + (Y(0) + Z(1) + X(0) * Z(2), "qwc"), + (Y(0) + Z(1) + X(0) * Z(2), "gc"), ], ) -def test_v_expectation_vmsa(hamiltonian, grouping): +def test_v_expectation_vmsa(terms, grouping): """Small scale tests of variance-based expectation value evaluation""" - n_qubits = 3 - circuit = Circuit(n_qubits) - circuit.add(gates.RX(_i, 0.1 * _i) for _i in range(n_qubits)) - circuit.add(gates.CNOT(_i, _i + 1) for _i in range(n_qubits - 1)) - circuit.add(gates.RZ(_i, 0.2 * _i) for _i in range(n_qubits)) + nqubits = 3 + circuit = Circuit(nqubits) + circuit.add(gates.RX(_i, 0.1 * _i) for _i in range(nqubits)) + circuit.add(gates.CNOT(_i, _i + 1) for _i in range(nqubits - 1)) + circuit.add(gates.RZ(_i, 0.2 * _i) for _i in range(nqubits)) + hamiltonian = SymbolicHamiltonian(terms, nqubits=nqubits) expected = hamiltonian.expectation(circuit) - n_shots = 50000 + n_shots = 50_000 n_trial_shots = 2000 test = v_expectation( circuit, From 6035d854064f4b02b48b586a32155fda99b2804a Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 3 Aug 2026 17:25:06 +0800 Subject: [PATCH 13/36] Minor tidying edits to code --- src/qibochem/measurement/optimization.py | 2 +- src/qibochem/measurement/util.py | 34 +++++++++++------------- 2 files changed, 16 insertions(+), 20 deletions(-) diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index c0ce366e..f021c879 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -135,7 +135,7 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl if dim_v != dim_symplectic: nullspace = _binary_nullspace(v_basis) # Interchange the 1st/2nd half of the indices to get nullspace in a symplectic sense - nullspace = np.concatenate((nullspace[:, dim_symplectic:], nullspace[:, :dim_symplectic]), axis=1) + nullspace = nullspace[:, np.r_[dim_symplectic : 2 * dim_symplectic, 0:dim_symplectic]] nullspace = _binary_gaussian_elimination(nullspace) v_basis = _lagrangian_subspace(nullspace) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index c81c93c5..9aa88994 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -5,7 +5,6 @@ import networkx as nx import numpy as np from qibo import gates -from qibo.config import raise_error # Mapping of Pauli operators to a symplectic (binary) representation, folowing the convention of (X|Z) PAULI_BINARY = {"I": (0, 0), "X": (1, 0), "Y": (1, 1), "Z": (0, 1)} @@ -167,28 +166,25 @@ def _binary_nullspace(binary_matrix: np.ndarray) -> np.ndarray: def _lagrangian_subspace(vector_space: np.ndarray) -> np.ndarray: """Find Lagrangian subspace of the given vector space; the symplectic nullspace in this context""" - cp_vector_space = np.array(vector_space) - # While loop to remove rows from cp_vector_space until cp_vector_space.shape matches (N, 2N) - while True: + # Remove rows from cp_vector_space until cp_vector_space.shape matches (N, 2N) + while vector_space.shape[0] > (vector_space.shape[1] // 2): anticommuting_vector_indices, anticommuting_vectors = None, None # Find a pair of anti-commuting vectors in vector_space - for _i1, _v1 in enumerate(cp_vector_space): - for _i2, _v2 in enumerate(cp_vector_space): - if _i2 > _i1 and _symplectic_inner_product(_v1, _v2) == 1: - anticommuting_vector_indices = [_i1, _i2] - anticommuting_vectors = cp_vector_space[anticommuting_vector_indices] + for i1, v1 in enumerate(vector_space): + for i2, v2 in enumerate(vector_space): + if i2 > i1 and _symplectic_inner_product(v1, v2) == 1: + anticommuting_vector_indices = [i1, i2] + anticommuting_vectors = vector_space[anticommuting_vector_indices] break if anticommuting_vector_indices is not None: break - if cp_vector_space.shape[0] == (cp_vector_space.shape[1] // 2): - break # Remove the two anti-commuting vectors from the basis - space_to_orthogonalize = np.delete(cp_vector_space, anticommuting_vector_indices, axis=0) - for _i1, vector in enumerate(space_to_orthogonalize): - for _i2, anticommuting_vector in enumerate(anticommuting_vectors): - space_to_orthogonalize[_i1] += ( - _symplectic_inner_product(vector, anticommuting_vectors[1 - _i2]) * anticommuting_vector + space_to_orthogonalize = np.delete(vector_space, anticommuting_vector_indices, axis=0) + for i1, vector in enumerate(space_to_orthogonalize): + for i2, anticommuting_vector in enumerate(anticommuting_vectors): + space_to_orthogonalize[i1] += ( + _symplectic_inner_product(vector, anticommuting_vectors[1 - i2]) * anticommuting_vector ) space_to_orthogonalize = space_to_orthogonalize % 2 @@ -196,10 +192,10 @@ def _lagrangian_subspace(vector_space: np.ndarray) -> np.ndarray: first_nonzero_col = np.argmax(anticommuting_vectors, axis=1) selected_vector = anticommuting_vectors[np.argmax(first_nonzero_col)] - cp_vector_space = np.append([selected_vector], space_to_orthogonalize, axis=0) - cp_vector_space = _binary_gaussian_elimination(cp_vector_space) + vector_space = np.append([selected_vector], space_to_orthogonalize, axis=0) + vector_space = _binary_gaussian_elimination(vector_space) - return cp_vector_space + return vector_space def _sort_tau_terms(v_basis: np.ndarray) -> np.ndarray: From 5f7214abb9d9ca2da9569383c691d3b36755495f Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Tue, 4 Aug 2026 14:23:55 +0800 Subject: [PATCH 14/36] Adding another test --- tests/test_expectation_samples.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_expectation_samples.py b/tests/test_expectation_samples.py index effe6495..6bfb98a6 100644 --- a/tests/test_expectation_samples.py +++ b/tests/test_expectation_samples.py @@ -106,6 +106,7 @@ def test_measurement_grouping_functionality(grouping, terms): (gates.RX(0, theta=-np.pi / 2), gates.H(1)), ), (Y(1) * Y(2) + X(0) * X(1) * Z(2), 3, (gates.H(0), gates.H(1))), + (Y(0) * X(1) + X(0) * Y(1) * Z(2), 3, (gates.H(0), gates.S(0), gates.H(1), gates.X(2), gates.H(2))), ], ) def test_measurement_grouping_extra_tests(grouping, terms, nqubits, gates_to_add): From 7b375f22bb319b11fa40592537f740b9422be524 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Wed, 5 Aug 2026 15:45:29 +0800 Subject: [PATCH 15/36] Still stuck on this... --- src/qibochem/measurement/optimization.py | 5 +++++ src/qibochem/measurement/util.py | 15 +++++++++------ 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index f021c879..6ef83403 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -139,15 +139,20 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl nullspace = _binary_gaussian_elimination(nullspace) v_basis = _lagrangian_subspace(nullspace) + # print("v_basis:\n", v_basis) + # print("Basis terms:", [_symplectic_to_pauli(x) for x in v_basis]) # Different methods of circuit synthesis if method == "chong": x_result = _solve_linear_system(v_basis, v_subspace) + # print(f"{x_result = }") # Map the solution onto the original set of qubits qubit_map = dict(zip(sorted({q for pauli_op in x_result for q in pauli_op}), term_qubits)) phase_factors = [_phase_factor(v_basis[pauli_op]) for pauli_op in x_result] + # print(f"{phase_factors = }") u_gates = _synthesise_circuit(v_basis) u_gates += [gates.SWAP(i, j) for i, j in qubit_map.items() if i != j] mapping = { + # term: phase * prod(Z(_i) for _i in soln) term: phase * prod(Z(qubit_map[_i]) for _i in soln) for term, phase, soln in zip(term_list, phase_factors, x_result) } diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index 9aa88994..c019aea2 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -304,15 +304,15 @@ def _make_x_matrix_full_rank(stabiliser_matrix: np.ndarray) -> list[gates.Gate]: z_matrix = stabiliser_matrix[:, dim_space:] # Need to find full rank submatrix in Z matrix for each of the zero rows in the X matrix + qubits = [] zero_row_indices = np.where(np.all(x_matrix == 0, axis=1))[0] - prev_swap = [] # To avoid swapping the same column twice while zero_row_indices.size > 0: # Select the first possible column for the first zero row for qubit in np.nonzero(z_matrix[zero_row_indices[0], :])[0]: - if qubit not in prev_swap: + if qubit not in qubits: stabiliser_matrix[:, [qubit, qubit + dim_space]] = stabiliser_matrix[:, [qubit + dim_space, qubit]] gates_list.append(gates.H(qubit)) - prev_swap.append(qubit) + qubits.append(qubit) break zero_row_indices = np.where(np.all(x_matrix == 0, axis=1))[0] return gates_list @@ -374,13 +374,13 @@ def _zero_z_matrix(stabiliser_matrix: np.ndarray) -> list[gates.Gate]: """ s_gates = [] cz_gates = [] - dim, _dim_space = stabiliser_matrix.shape - dim_space = _dim_space // 2 + dim, dim_space = stabiliser_matrix.shape + dim_space = dim_space // 2 # Following the algorithm in the paper, zero out the diagonal entries first for _i in range(dim): if stabiliser_matrix[_i, dim_space + _i] == 1: stabiliser_matrix[_i, dim_space + _i] = 0 - s_gates.append(gates.S(_i).dagger()) # Paper says S gate, but should be S.dagger? + s_gates.append(gates.S(_i)) # Then remove the off-diagonal terms in each row for _j in range(dim_space): if _j > _i and stabiliser_matrix[_i, dim_space + _j] == 1: @@ -403,10 +403,13 @@ def _synthesise_circuit(v_basis: np.ndarray) -> list[gates.Gate]: rotation_gates = [] # 1. Apply H gates to transform 'X matrix' to full rank rotation_gates += _make_x_matrix_full_rank(stabiliser_matrix) + # print("Matrix:\n", stabiliser_matrix) # 2. Row-reduce 'X matrix' to I using CNOT/SWAP gates rotation_gates += _col_reduce_x_matrix(stabiliser_matrix) + # print("Matrix:\n", stabiliser_matrix) # 3. Remove all non-zero entries on 'Z matrix' using S and CZ gates rotation_gates += _zero_z_matrix(stabiliser_matrix) + # print("Matrix:\n", stabiliser_matrix) # 4. Apply H to each qubit to swap the 'X' and 'Z' matrices rotation_gates += [gates.H(_i) for _i in range(n_qubits)] return rotation_gates From 057f7ffd9a14a7c56a30b12f678788331d0c9ee1 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Wed, 5 Aug 2026 15:51:01 +0800 Subject: [PATCH 16/36] Fix single test --- tests/test_measurement_util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_measurement_util.py b/tests/test_measurement_util.py index 5763c31a..7cff4d2e 100644 --- a/tests/test_measurement_util.py +++ b/tests/test_measurement_util.py @@ -274,4 +274,4 @@ def test_zero_z_matrix(): ) gates_list = _zero_z_matrix(stabiliser_matrix) # Single column operation, should have only CNOT gate - assert len(gates_list) == 1 and gates_list[0].name == "sdg" + assert len(gates_list) == 1 and gates_list[0].name == "s" From 18381eb5142aeb1c0fc6aeb05a3496cf0842b9c6 Mon Sep 17 00:00:00 2001 From: shangtai Date: Sun, 9 Aug 2026 01:51:26 +0800 Subject: [PATCH 17/36] attempt to fix the phase --- src/qibochem/measurement/util.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index c019aea2..f9375f01 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -260,14 +260,14 @@ def _solve_linear_system(binary_matrix: np.ndarray, vector: np.ndarray) -> list[ def _single_qubit_phase_factor(pauli_ops: list[np.ndarray]) -> complex: """Compute the phase factor w.r.t. the product of multiple Pauli operators for a single qubit""" # Initialise as 1.0*I, then multiply with each Pauli operator acting on that qubit - coeff, current_pauli_op = 1.0, np.zeros(2) + coeff, current_pauli_op = 1.0, None for pauli_op in pauli_ops: # If I, just skip - if SYMPLECTIC_INDEX[tuple(current_pauli_op)] == 0: - current_pauli_op = pauli_op - continue if SYMPLECTIC_INDEX[tuple(pauli_op)] == 0: continue + if current_pauli_op is None: + current_pauli_op = pauli_op + continue # Multiply by some phase factor depending on what Pauli operators are involved coeff *= SYMPLECTIC_PHASE_TABLE[SYMPLECTIC_INDEX[tuple(pauli_op)] - SYMPLECTIC_INDEX[tuple(current_pauli_op)]] current_pauli_op = (current_pauli_op + pauli_op) % 2 @@ -276,10 +276,6 @@ def _single_qubit_phase_factor(pauli_ops: list[np.ndarray]) -> complex: def _phase_factor(pauli_terms: list[np.ndarray]) -> int: """Compute phase factor of a product of mutually commuting Pauli terms (in symplectic form). Returns: 1 or -1""" - # Singleton case is trivial: 1 - if len(pauli_terms) == 1: - return 1 - # >1 term: dim = pauli_terms[0].shape[0] // 2 coefficient = 1.0 for qubit in range(dim): From 1b64f8a7d9cf4586de078f5b38dc7ae1de86592b Mon Sep 17 00:00:00 2001 From: shangtai Date: Sun, 9 Aug 2026 02:18:28 +0800 Subject: [PATCH 18/36] attempt to fix the phase --- src/qibochem/measurement/util.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index f9375f01..053b1403 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -256,20 +256,27 @@ def _solve_linear_system(binary_matrix: np.ndarray, vector: np.ndarray) -> list[ # Get non-zero entries in each column on RHS of rref_aug_matrix => Solution for respective vector in b return [np.nonzero(rref_aug_matrix[:, binary_matrix.shape[0] + i])[0].tolist() for i in range(vector.shape[0])] +PAULI_MULTIPLICATION_PHASE = { + ((1, 0), (0, 1)): 1.0j, # X * Z = iY + ((0, 1), (1, 0)): -1.0j, # Z * X = -iY + ((1, 0), (1, 1)): -1.0j, # X * Y = -iZ + ((1, 1), (1, 0)): 1.0j, # Y * X = iZ + ((0, 1), (1, 1)): 1.0j, # Z * Y = iX + ((1, 1), (0, 1)): -1.0j, # Y * Z = -iX +} def _single_qubit_phase_factor(pauli_ops: list[np.ndarray]) -> complex: - """Compute the phase factor w.r.t. the product of multiple Pauli operators for a single qubit""" - # Initialise as 1.0*I, then multiply with each Pauli operator acting on that qubit coeff, current_pauli_op = 1.0, None for pauli_op in pauli_ops: - # If I, just skip - if SYMPLECTIC_INDEX[tuple(pauli_op)] == 0: + key = tuple(pauli_op) + if SYMPLECTIC_INDEX[key] == 0: continue if current_pauli_op is None: current_pauli_op = pauli_op continue - # Multiply by some phase factor depending on what Pauli operators are involved - coeff *= SYMPLECTIC_PHASE_TABLE[SYMPLECTIC_INDEX[tuple(pauli_op)] - SYMPLECTIC_INDEX[tuple(current_pauli_op)]] + + current_key = tuple(current_pauli_op) + coeff *= PAULI_MULTIPLICATION_PHASE[(current_key, key)] current_pauli_op = (current_pauli_op + pauli_op) % 2 return coeff From 107250c4f5a1ecc1f68c9c8ec0b801be9d694538 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 8 Aug 2026 18:18:44 +0000 Subject: [PATCH 19/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/qibochem/measurement/util.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index 053b1403..7e304825 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -256,15 +256,17 @@ def _solve_linear_system(binary_matrix: np.ndarray, vector: np.ndarray) -> list[ # Get non-zero entries in each column on RHS of rref_aug_matrix => Solution for respective vector in b return [np.nonzero(rref_aug_matrix[:, binary_matrix.shape[0] + i])[0].tolist() for i in range(vector.shape[0])] + PAULI_MULTIPLICATION_PHASE = { - ((1, 0), (0, 1)): 1.0j, # X * Z = iY + ((1, 0), (0, 1)): 1.0j, # X * Z = iY ((0, 1), (1, 0)): -1.0j, # Z * X = -iY ((1, 0), (1, 1)): -1.0j, # X * Y = -iZ - ((1, 1), (1, 0)): 1.0j, # Y * X = iZ - ((0, 1), (1, 1)): 1.0j, # Z * Y = iX + ((1, 1), (1, 0)): 1.0j, # Y * X = iZ + ((0, 1), (1, 1)): 1.0j, # Z * Y = iX ((1, 1), (0, 1)): -1.0j, # Y * Z = -iX } + def _single_qubit_phase_factor(pauli_ops: list[np.ndarray]) -> complex: coeff, current_pauli_op = 1.0, None for pauli_op in pauli_ops: From a0875a58d4ad5300adc8163feca563d9d7d58f5d Mon Sep 17 00:00:00 2001 From: shangtai Date: Sun, 9 Aug 2026 02:23:52 +0800 Subject: [PATCH 20/36] attempt to fix the phase --- src/qibochem/measurement/util.py | 31 +++++++++++++------------------ 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index 7e304825..c019aea2 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -257,34 +257,29 @@ def _solve_linear_system(binary_matrix: np.ndarray, vector: np.ndarray) -> list[ return [np.nonzero(rref_aug_matrix[:, binary_matrix.shape[0] + i])[0].tolist() for i in range(vector.shape[0])] -PAULI_MULTIPLICATION_PHASE = { - ((1, 0), (0, 1)): 1.0j, # X * Z = iY - ((0, 1), (1, 0)): -1.0j, # Z * X = -iY - ((1, 0), (1, 1)): -1.0j, # X * Y = -iZ - ((1, 1), (1, 0)): 1.0j, # Y * X = iZ - ((0, 1), (1, 1)): 1.0j, # Z * Y = iX - ((1, 1), (0, 1)): -1.0j, # Y * Z = -iX -} - - def _single_qubit_phase_factor(pauli_ops: list[np.ndarray]) -> complex: - coeff, current_pauli_op = 1.0, None + """Compute the phase factor w.r.t. the product of multiple Pauli operators for a single qubit""" + # Initialise as 1.0*I, then multiply with each Pauli operator acting on that qubit + coeff, current_pauli_op = 1.0, np.zeros(2) for pauli_op in pauli_ops: - key = tuple(pauli_op) - if SYMPLECTIC_INDEX[key] == 0: - continue - if current_pauli_op is None: + # If I, just skip + if SYMPLECTIC_INDEX[tuple(current_pauli_op)] == 0: current_pauli_op = pauli_op continue - - current_key = tuple(current_pauli_op) - coeff *= PAULI_MULTIPLICATION_PHASE[(current_key, key)] + if SYMPLECTIC_INDEX[tuple(pauli_op)] == 0: + continue + # Multiply by some phase factor depending on what Pauli operators are involved + coeff *= SYMPLECTIC_PHASE_TABLE[SYMPLECTIC_INDEX[tuple(pauli_op)] - SYMPLECTIC_INDEX[tuple(current_pauli_op)]] current_pauli_op = (current_pauli_op + pauli_op) % 2 return coeff def _phase_factor(pauli_terms: list[np.ndarray]) -> int: """Compute phase factor of a product of mutually commuting Pauli terms (in symplectic form). Returns: 1 or -1""" + # Singleton case is trivial: 1 + if len(pauli_terms) == 1: + return 1 + # >1 term: dim = pauli_terms[0].shape[0] // 2 coefficient = 1.0 for qubit in range(dim): From 89d394ad86a66f205914ff504f53e5ed58b0502e Mon Sep 17 00:00:00 2001 From: shangtai Date: Sun, 9 Aug 2026 02:57:12 +0800 Subject: [PATCH 21/36] attempt to fix the phase --- src/qibochem/measurement/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index c019aea2..4a783647 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -10,7 +10,7 @@ PAULI_BINARY = {"I": (0, 0), "X": (1, 0), "Y": (1, 1), "Z": (0, 1)} BINARY_PAULI = {symplectic: pauli for pauli, symplectic in PAULI_BINARY.items()} -SYMPLECTIC_PHASE_TABLE = [1.0, 1.0j, -1.0j] +SYMPLECTIC_PHASE_TABLE = [1.0, -1.0j, 1.0j] SYMPLECTIC_INDEX = {symplectic: index for index, symplectic in enumerate(BINARY_PAULI.keys())} From a803385232cded8f2875d2004033002a16470790 Mon Sep 17 00:00:00 2001 From: shangtai Date: Sun, 9 Aug 2026 03:19:04 +0800 Subject: [PATCH 22/36] attempt to fix the phase --- src/qibochem/measurement/util.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index 4a783647..b6ca63c4 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -260,13 +260,16 @@ def _solve_linear_system(binary_matrix: np.ndarray, vector: np.ndarray) -> list[ def _single_qubit_phase_factor(pauli_ops: list[np.ndarray]) -> complex: """Compute the phase factor w.r.t. the product of multiple Pauli operators for a single qubit""" # Initialise as 1.0*I, then multiply with each Pauli operator acting on that qubit - coeff, current_pauli_op = 1.0, np.zeros(2) + coeff, current_pauli_op = 1.0, None for pauli_op in pauli_ops: # If I, just skip - if SYMPLECTIC_INDEX[tuple(current_pauli_op)] == 0: + if SYMPLECTIC_INDEX[tuple(pauli_op)] == 0: + continue + if current_pauli_op is None: current_pauli_op = pauli_op continue - if SYMPLECTIC_INDEX[tuple(pauli_op)] == 0: + if SYMPLECTIC_INDEX[tuple(current_pauli_op)] == 0: + current_pauli_op = pauli_op continue # Multiply by some phase factor depending on what Pauli operators are involved coeff *= SYMPLECTIC_PHASE_TABLE[SYMPLECTIC_INDEX[tuple(pauli_op)] - SYMPLECTIC_INDEX[tuple(current_pauli_op)]] From 224d65a3707f1d7f095d0efac84687d4e111fb57 Mon Sep 17 00:00:00 2001 From: shangtai Date: Mon, 10 Aug 2026 03:42:13 +0800 Subject: [PATCH 23/36] attempt to fix the phase --- src/qibochem/measurement/util.py | 28 +++++++--------------------- 1 file changed, 7 insertions(+), 21 deletions(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index b6ca63c4..c8d83810 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -10,7 +10,7 @@ PAULI_BINARY = {"I": (0, 0), "X": (1, 0), "Y": (1, 1), "Z": (0, 1)} BINARY_PAULI = {symplectic: pauli for pauli, symplectic in PAULI_BINARY.items()} -SYMPLECTIC_PHASE_TABLE = [1.0, -1.0j, 1.0j] +SYMPLECTIC_PHASE_TABLE = [1.0, 1.0j, -1.0j] SYMPLECTIC_INDEX = {symplectic: index for index, symplectic in enumerate(BINARY_PAULI.keys())} @@ -259,34 +259,20 @@ def _solve_linear_system(binary_matrix: np.ndarray, vector: np.ndarray) -> list[ def _single_qubit_phase_factor(pauli_ops: list[np.ndarray]) -> complex: """Compute the phase factor w.r.t. the product of multiple Pauli operators for a single qubit""" - # Initialise as 1.0*I, then multiply with each Pauli operator acting on that qubit - coeff, current_pauli_op = 1.0, None + exponent = 0 # power of i, tracked mod 4 + r_x, r_z = 0, 0 # running product in symplectic form, starts at I = (0, 0) for pauli_op in pauli_ops: - # If I, just skip - if SYMPLECTIC_INDEX[tuple(pauli_op)] == 0: - continue - if current_pauli_op is None: - current_pauli_op = pauli_op - continue - if SYMPLECTIC_INDEX[tuple(current_pauli_op)] == 0: - current_pauli_op = pauli_op - continue - # Multiply by some phase factor depending on what Pauli operators are involved - coeff *= SYMPLECTIC_PHASE_TABLE[SYMPLECTIC_INDEX[tuple(pauli_op)] - SYMPLECTIC_INDEX[tuple(current_pauli_op)]] - current_pauli_op = (current_pauli_op + pauli_op) % 2 - return coeff + b_x, b_z = int(pauli_op[0]), int(pauli_op[1]) + exponent = (exponent + r_x * b_z - r_z * b_x) % 4 + r_x, r_z = r_x ^ b_x, r_z ^ b_z + return 1j ** exponent def _phase_factor(pauli_terms: list[np.ndarray]) -> int: """Compute phase factor of a product of mutually commuting Pauli terms (in symplectic form). Returns: 1 or -1""" - # Singleton case is trivial: 1 - if len(pauli_terms) == 1: - return 1 - # >1 term: dim = pauli_terms[0].shape[0] // 2 coefficient = 1.0 for qubit in range(dim): - # Get all Pauli operators for a particular qubit pauli_ops = [pauli_term[[qubit, qubit + dim]] for pauli_term in pauli_terms] coefficient *= _single_qubit_phase_factor(pauli_ops) return int(np.real_if_close(coefficient)) From 9e52c1adaaa8004208d69d9ae945d0e830a32a7e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 9 Aug 2026 19:42:58 +0000 Subject: [PATCH 24/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/qibochem/measurement/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index c8d83810..1add7a08 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -265,7 +265,7 @@ def _single_qubit_phase_factor(pauli_ops: list[np.ndarray]) -> complex: b_x, b_z = int(pauli_op[0]), int(pauli_op[1]) exponent = (exponent + r_x * b_z - r_z * b_x) % 4 r_x, r_z = r_x ^ b_x, r_z ^ b_z - return 1j ** exponent + return 1j**exponent def _phase_factor(pauli_terms: list[np.ndarray]) -> int: From d801295ff50e2ac1e07899d4ee25fee3b243414d Mon Sep 17 00:00:00 2001 From: shangtai Date: Mon, 10 Aug 2026 03:53:16 +0800 Subject: [PATCH 25/36] attempt to fix the phase --- src/qibochem/measurement/util.py | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index c8d83810..c019aea2 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -259,20 +259,31 @@ def _solve_linear_system(binary_matrix: np.ndarray, vector: np.ndarray) -> list[ def _single_qubit_phase_factor(pauli_ops: list[np.ndarray]) -> complex: """Compute the phase factor w.r.t. the product of multiple Pauli operators for a single qubit""" - exponent = 0 # power of i, tracked mod 4 - r_x, r_z = 0, 0 # running product in symplectic form, starts at I = (0, 0) + # Initialise as 1.0*I, then multiply with each Pauli operator acting on that qubit + coeff, current_pauli_op = 1.0, np.zeros(2) for pauli_op in pauli_ops: - b_x, b_z = int(pauli_op[0]), int(pauli_op[1]) - exponent = (exponent + r_x * b_z - r_z * b_x) % 4 - r_x, r_z = r_x ^ b_x, r_z ^ b_z - return 1j ** exponent + # If I, just skip + if SYMPLECTIC_INDEX[tuple(current_pauli_op)] == 0: + current_pauli_op = pauli_op + continue + if SYMPLECTIC_INDEX[tuple(pauli_op)] == 0: + continue + # Multiply by some phase factor depending on what Pauli operators are involved + coeff *= SYMPLECTIC_PHASE_TABLE[SYMPLECTIC_INDEX[tuple(pauli_op)] - SYMPLECTIC_INDEX[tuple(current_pauli_op)]] + current_pauli_op = (current_pauli_op + pauli_op) % 2 + return coeff def _phase_factor(pauli_terms: list[np.ndarray]) -> int: """Compute phase factor of a product of mutually commuting Pauli terms (in symplectic form). Returns: 1 or -1""" + # Singleton case is trivial: 1 + if len(pauli_terms) == 1: + return 1 + # >1 term: dim = pauli_terms[0].shape[0] // 2 coefficient = 1.0 for qubit in range(dim): + # Get all Pauli operators for a particular qubit pauli_ops = [pauli_term[[qubit, qubit + dim]] for pauli_term in pauli_terms] coefficient *= _single_qubit_phase_factor(pauli_ops) return int(np.real_if_close(coefficient)) From 96e82b48634441f00a2d125634b058983d96ea16 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Tue, 18 Aug 2026 11:26:46 +0800 Subject: [PATCH 26/36] Minor cleaning of test_molecule --- tests/test_molecule.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/test_molecule.py b/tests/test_molecule.py index 4a028670..54132833 100644 --- a/tests/test_molecule.py +++ b/tests/test_molecule.py @@ -7,7 +7,6 @@ import numpy as np import openfermion import pytest -from qibo import gates, models from qibo.hamiltonians import SymbolicHamiltonian from qibo.symbols import X, Z @@ -29,7 +28,7 @@ def test_pyscf_driver(xyz_file, expected): file_path = Path("tests", "data") / Path(xyz_file) # In case .xyz files somehow not found if not file_path.is_file(): - with open(file_path, "w") as file_handler: + with open(file_path, "w", encoding="utf-8") as file_handler: if xyz_file == "lih.xyz": file_handler.write("2\n0 1\nLi 0.0 0.0 0.0\nH 0.0 0.0 1.2\n") elif xyz_file == "h2.xyz": From 7b8f226494ac4d03a4bc09d1006f5330c6f480aa Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Tue, 18 Aug 2026 11:32:27 +0800 Subject: [PATCH 27/36] Add array to track phase of basis terms --- src/qibochem/measurement/optimization.py | 11 +--- src/qibochem/measurement/util.py | 71 ++++++++++++++---------- 2 files changed, 44 insertions(+), 38 deletions(-) diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index 6ef83403..07264206 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -139,22 +139,17 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl nullspace = _binary_gaussian_elimination(nullspace) v_basis = _lagrangian_subspace(nullspace) - # print("v_basis:\n", v_basis) - # print("Basis terms:", [_symplectic_to_pauli(x) for x in v_basis]) # Different methods of circuit synthesis if method == "chong": x_result = _solve_linear_system(v_basis, v_subspace) - # print(f"{x_result = }") # Map the solution onto the original set of qubits qubit_map = dict(zip(sorted({q for pauli_op in x_result for q in pauli_op}), term_qubits)) phase_factors = [_phase_factor(v_basis[pauli_op]) for pauli_op in x_result] - # print(f"{phase_factors = }") - u_gates = _synthesise_circuit(v_basis) + u_gates, phases = _synthesise_circuit(v_basis) u_gates += [gates.SWAP(i, j) for i, j in qubit_map.items() if i != j] mapping = { - # term: phase * prod(Z(_i) for _i in soln) - term: phase * prod(Z(qubit_map[_i]) for _i in soln) - for term, phase, soln in zip(term_list, phase_factors, x_result) + term: circuit_phase * phase * prod(Z(qubit_map[_i]) for _i in soln) + for term, circuit_phase, phase, soln in zip(term_list, phases, phase_factors, x_result) } elif method == "izmaylov": v_basis = _sort_tau_terms(v_basis) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index c019aea2..b99b9664 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -289,9 +289,9 @@ def _phase_factor(pauli_terms: list[np.ndarray]) -> int: return int(np.real_if_close(coefficient)) -def _make_x_matrix_full_rank(stabiliser_matrix: np.ndarray) -> list[gates.Gate]: +def _make_x_matrix_full_rank(stabiliser_matrix: np.ndarray, phases: np.ndarray) -> list[gates.Gate]: """ - Modifies stabiliser_matrix (in-place) to transform 'X matrix' to full rank, with H gates representing each 'swap' + Modifies stabiliser_matrix and phases in-place to transform 'X matrix' to full rank, with H gates representing each 'swap' of columns between the 'Z' and 'X' matrices. Note: stabiliser_matrix should already be in reduced row echelon form Returns: @@ -310,6 +310,8 @@ def _make_x_matrix_full_rank(stabiliser_matrix: np.ndarray) -> list[gates.Gate]: # Select the first possible column for the first zero row for qubit in np.nonzero(z_matrix[zero_row_indices[0], :])[0]: if qubit not in qubits: + # For S(a)/H(a): r_i := r_i + x_{i,a} z_{i,a} for all i + phases += stabiliser_matrix[:, qubit] * stabiliser_matrix[:, qubit + dim_space] stabiliser_matrix[:, [qubit, qubit + dim_space]] = stabiliser_matrix[:, [qubit + dim_space, qubit]] gates_list.append(gates.H(qubit)) qubits.append(qubit) @@ -318,9 +320,9 @@ def _make_x_matrix_full_rank(stabiliser_matrix: np.ndarray) -> list[gates.Gate]: return gates_list -def _col_reduce_x_matrix(stabiliser_matrix: np.ndarray) -> list[gates.Gate]: +def _col_reduce_x_matrix(stabiliser_matrix: np.ndarray, phases: np.ndarray) -> list[gates.Gate]: """ - Modifies stabiliser_matrix in-place to transform the X matrix to I, using CNOT/SWAP gates + Modifies stabiliser_matrix and phases in-place to transform the X matrix to I, using CNOT/SWAP gates Returns: list[gates.Gate]: List of CNOT/SWAP gates to be added to the circuit @@ -352,22 +354,27 @@ def _col_reduce_x_matrix(stabiliser_matrix: np.ndarray) -> list[gates.Gate]: # Remove all nonzero entries on row _i using CNOT gates for col in nonzero_cols: + # For CNOT(a, b): r_i := r_i + x_{i,a} z_{i,b} (x_{i,b} + z_{i,a} + 1), for all i + phase_changes = stabiliser_matrix[:, col] + stabiliser_matrix[:, pivot_col] + 1 + phase_changes %= 2 + phase_changes *= stabiliser_matrix[:, pivot_col] * stabiliser_matrix[:, col] + phases += phase_changes # X matrix: Add pivot column to column with 1 stabiliser_matrix[:, col] += stabiliser_matrix[:, pivot_col] # Z matrix: Add (column with 1)^th column to pivot column stabiliser_matrix[:, pivot_col + dim_space] += stabiliser_matrix[:, col + dim_space] stabiliser_matrix %= 2 - gates_list.append(gates.CNOT(col, pivot_col)) + gates_list.append(gates.CNOT(pivot_col, col)) pivot_col += 1 return gates_list -def _zero_z_matrix(stabiliser_matrix: np.ndarray) -> list[gates.Gate]: +def _zero_z_matrix(stabiliser_matrix: np.ndarray, phases: np.ndarray) -> list[gates.Gate]: """ - Modifies stabiliser_matrix in-place to transform the Z matrix to a zero matrix. - 1. S gates used to set diagonal entries on Z matrix - 2. CZ gates used to remove off-diagonal entries on Z matrix + Modifies stabiliser_matrix and phases in-place to transform the Z matrix to a zero matrix. + 1. S gates used to set diagonal entries on Z matrix (Phases updated) + 2. CZ gates used to remove off-diagonal entries on Z matrix (Phases not updated) Returns: list[gates.Gate]: List of S and CZ gates to be added to the circuit @@ -377,39 +384,43 @@ def _zero_z_matrix(stabiliser_matrix: np.ndarray) -> list[gates.Gate]: dim, dim_space = stabiliser_matrix.shape dim_space = dim_space // 2 # Following the algorithm in the paper, zero out the diagonal entries first - for _i in range(dim): - if stabiliser_matrix[_i, dim_space + _i] == 1: - stabiliser_matrix[_i, dim_space + _i] = 0 - s_gates.append(gates.S(_i)) + for i in range(dim): + if stabiliser_matrix[i, dim_space + i] == 1: + # For S(a)/H(a): r_i := r_i + x_{i,a} z_{i,a} for all i + phases += stabiliser_matrix[:, i] * stabiliser_matrix[:, i + dim_space] + stabiliser_matrix[i, dim_space + i] = 0 + s_gates.append(gates.S(i)) # Then remove the off-diagonal terms in each row - for _j in range(dim_space): - if _j > _i and stabiliser_matrix[_i, dim_space + _j] == 1: - stabiliser_matrix[_i, dim_space + _j] = 0 - stabiliser_matrix[_j, dim_space + _i] = 0 - cz_gates.append(gates.CZ(_i, _j)) + for j in range(dim_space): + if j > i and stabiliser_matrix[i, dim_space + j] == 1: + # Note: Not updating the phases w.r.t. CZ + stabiliser_matrix[i, dim_space + j] = 0 + stabiliser_matrix[j, dim_space + i] = 0 + cz_gates.append(gates.CZ(i, j)) return s_gates + cz_gates -def _synthesise_circuit(v_basis: np.ndarray) -> list[gates.Gate]: +def _synthesise_circuit(v_basis: np.ndarray) -> tuple[list[gates.Gate], list[int]]: """ Gets the basis rotation gates for rotating the initial measurement basis into the computational basis. The stabiliser matrix (v_basis) follows the format of (X|Z) matrices. Returns: list[gates.Gate]: Gates to be added after the circuit ansatz + list[int]: Phases of the measured basis terms """ - stabiliser_matrix = np.array(v_basis) - n_qubits = stabiliser_matrix.shape[1] // 2 + stabiliser_matrix = np.array(v_basis, dtype=np.uint8) + nqubits = stabiliser_matrix.shape[0] + phases = np.array([[0 for _ in range(nqubits)]], dtype=np.uint8) # To keep track of phases rotation_gates = [] # 1. Apply H gates to transform 'X matrix' to full rank - rotation_gates += _make_x_matrix_full_rank(stabiliser_matrix) - # print("Matrix:\n", stabiliser_matrix) + rotation_gates += _make_x_matrix_full_rank(stabiliser_matrix, phases) # 2. Row-reduce 'X matrix' to I using CNOT/SWAP gates - rotation_gates += _col_reduce_x_matrix(stabiliser_matrix) - # print("Matrix:\n", stabiliser_matrix) + rotation_gates += _col_reduce_x_matrix(stabiliser_matrix, phases) # 3. Remove all non-zero entries on 'Z matrix' using S and CZ gates - rotation_gates += _zero_z_matrix(stabiliser_matrix) - # print("Matrix:\n", stabiliser_matrix) - # 4. Apply H to each qubit to swap the 'X' and 'Z' matrices - rotation_gates += [gates.H(_i) for _i in range(n_qubits)] - return rotation_gates + rotation_gates += _zero_z_matrix(stabiliser_matrix, phases) + # 4. Apply H to each qubit to swap the 'X' and 'Z' matrices. Note: Not gonna update phases here + rotation_gates += [gates.H(i) for i in range(nqubits)] + # Update circuit phase factors to be 1 or -1 + phases = [-1 if x else 1 for x in phases[0]] + return rotation_gates, phases From b58da1154bb30d3319e63e1f7a94dd1c22e93fec Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Tue, 18 Aug 2026 11:32:58 +0800 Subject: [PATCH 28/36] Update test functions --- tests/test_measurement_util.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/tests/test_measurement_util.py b/tests/test_measurement_util.py index 7cff4d2e..5d268409 100644 --- a/tests/test_measurement_util.py +++ b/tests/test_measurement_util.py @@ -257,9 +257,11 @@ def test_col_reduce_x_matrix(): ], dtype=np.uint8, ) - gates_list = _col_reduce_x_matrix(stabiliser_matrix) + phases = np.array([0, 0], dtype=np.uint8) + gates_list = _col_reduce_x_matrix(stabiliser_matrix, phases) # Single column operation, should have only CNOT gate assert len(gates_list) == 1 and gates_list[0].name == "cx" + assert phases[0] == 1 def test_zero_z_matrix(): @@ -272,6 +274,9 @@ def test_zero_z_matrix(): ], dtype=np.uint8, ) - gates_list = _zero_z_matrix(stabiliser_matrix) + phases = np.array([0, 0, 0, 0], dtype=np.uint8) + gates_list = _zero_z_matrix(stabiliser_matrix, phases) + print(phases) # Single column operation, should have only CNOT gate assert len(gates_list) == 1 and gates_list[0].name == "s" + assert phases[0] == 1 From 554825208682afee12687f7f13dc49f034b15af7 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Tue, 18 Aug 2026 12:03:34 +0800 Subject: [PATCH 29/36] Fix code coverage --- tests/test_measurement_util.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/tests/test_measurement_util.py b/tests/test_measurement_util.py index 5d268409..78af1b61 100644 --- a/tests/test_measurement_util.py +++ b/tests/test_measurement_util.py @@ -262,6 +262,20 @@ def test_col_reduce_x_matrix(): # Single column operation, should have only CNOT gate assert len(gates_list) == 1 and gates_list[0].name == "cx" assert phases[0] == 1 + # Code coverage for Gaussian elimination. Note: Input matrix isn't a commuting set, so shouldn't ever need + control = np.array( + [ + [1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + ], + dtype=np.uint8, + ) + stabiliser_matrix = np.array(control, dtype=np.uint8) + phases = np.array([0, 0], dtype=np.uint8) + _gates = _col_reduce_x_matrix(stabiliser_matrix, phases) + # No change to stabiliser matrix + assert np.array_equal(control, stabiliser_matrix) def test_zero_z_matrix(): From 63f5b35626925e6552534771b50e81985603c400 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Wed, 19 Aug 2026 10:25:40 +0800 Subject: [PATCH 30/36] Bugfix: Phase tracking - Product of phases was wrong - Wrong column reference for CNOT gate phase update --- src/qibochem/measurement/optimization.py | 4 ++-- src/qibochem/measurement/util.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index 07264206..88ba0725 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -148,8 +148,8 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl u_gates, phases = _synthesise_circuit(v_basis) u_gates += [gates.SWAP(i, j) for i, j in qubit_map.items() if i != j] mapping = { - term: circuit_phase * phase * prod(Z(qubit_map[_i]) for _i in soln) - for term, circuit_phase, phase, soln in zip(term_list, phases, phase_factors, x_result) + term: phase * prod(phases[i] * Z(qubit_map[i]) for i in soln) + for term, phase, soln in zip(term_list, phase_factors, x_result) } elif method == "izmaylov": v_basis = _sort_tau_terms(v_basis) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index b99b9664..ff00760b 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -355,9 +355,9 @@ def _col_reduce_x_matrix(stabiliser_matrix: np.ndarray, phases: np.ndarray) -> l # Remove all nonzero entries on row _i using CNOT gates for col in nonzero_cols: # For CNOT(a, b): r_i := r_i + x_{i,a} z_{i,b} (x_{i,b} + z_{i,a} + 1), for all i - phase_changes = stabiliser_matrix[:, col] + stabiliser_matrix[:, pivot_col] + 1 + phase_changes = stabiliser_matrix[:, col] + stabiliser_matrix[:, pivot_col + dim_space] + 1 phase_changes %= 2 - phase_changes *= stabiliser_matrix[:, pivot_col] * stabiliser_matrix[:, col] + phase_changes *= stabiliser_matrix[:, pivot_col] * stabiliser_matrix[:, col + dim_space] phases += phase_changes # X matrix: Add pivot column to column with 1 stabiliser_matrix[:, col] += stabiliser_matrix[:, pivot_col] @@ -373,7 +373,7 @@ def _col_reduce_x_matrix(stabiliser_matrix: np.ndarray, phases: np.ndarray) -> l def _zero_z_matrix(stabiliser_matrix: np.ndarray, phases: np.ndarray) -> list[gates.Gate]: """ Modifies stabiliser_matrix and phases in-place to transform the Z matrix to a zero matrix. - 1. S gates used to set diagonal entries on Z matrix (Phases updated) + 1. S gates used to set diagonal entries on Z matrix 2. CZ gates used to remove off-diagonal entries on Z matrix (Phases not updated) Returns: @@ -393,7 +393,7 @@ def _zero_z_matrix(stabiliser_matrix: np.ndarray, phases: np.ndarray) -> list[ga # Then remove the off-diagonal terms in each row for j in range(dim_space): if j > i and stabiliser_matrix[i, dim_space + j] == 1: - # Note: Not updating the phases w.r.t. CZ + # Note: Not updating phases here stabiliser_matrix[i, dim_space + j] = 0 stabiliser_matrix[j, dim_space + i] = 0 cz_gates.append(gates.CZ(i, j)) From 1716a482ed80f9914f8358e10df6d64ce6f6f93a Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Wed, 19 Aug 2026 10:32:24 +0800 Subject: [PATCH 31/36] Fix tests --- tests/test_expectation_samples.py | 1 + tests/test_measurement_util.py | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/test_expectation_samples.py b/tests/test_expectation_samples.py index 6bfb98a6..dc0f520a 100644 --- a/tests/test_expectation_samples.py +++ b/tests/test_expectation_samples.py @@ -99,6 +99,7 @@ def test_measurement_grouping_functionality(grouping, terms): @pytest.mark.parametrize( "terms,nqubits,gates_to_add", [ + (X(0) * X(1) + Y(0) * Z(1), 2, (gates.H(0), gates.H(1))), (0.5 * X(0) * Y(1) + Z(0) * Z(1) * Z(2), 3, (gates.H(0), gates.RX(1, theta=-np.pi / 2))), ( 0.5 * Y(0) * X(1) * Z(3) * Z(4) + Z(0) * Z(1) * X(2) * Z(3) * Z(4), diff --git a/tests/test_measurement_util.py b/tests/test_measurement_util.py index 78af1b61..fe743a72 100644 --- a/tests/test_measurement_util.py +++ b/tests/test_measurement_util.py @@ -259,9 +259,10 @@ def test_col_reduce_x_matrix(): ) phases = np.array([0, 0], dtype=np.uint8) gates_list = _col_reduce_x_matrix(stabiliser_matrix, phases) + print(phases) # Single column operation, should have only CNOT gate assert len(gates_list) == 1 and gates_list[0].name == "cx" - assert phases[0] == 1 + assert np.array_equal(phases, np.array([0, 0], dtype=np.uint8)) # Code coverage for Gaussian elimination. Note: Input matrix isn't a commuting set, so shouldn't ever need control = np.array( [ From 2e3773cc8a59058df48df5839efec088f6ba227b Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 31 Aug 2026 13:40:13 +0800 Subject: [PATCH 32/36] Use bitwise operations for phase tracking --- src/qibochem/measurement/util.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index ff00760b..34384ff8 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -311,7 +311,7 @@ def _make_x_matrix_full_rank(stabiliser_matrix: np.ndarray, phases: np.ndarray) for qubit in np.nonzero(z_matrix[zero_row_indices[0], :])[0]: if qubit not in qubits: # For S(a)/H(a): r_i := r_i + x_{i,a} z_{i,a} for all i - phases += stabiliser_matrix[:, qubit] * stabiliser_matrix[:, qubit + dim_space] + phases ^= stabiliser_matrix[:, qubit] * stabiliser_matrix[:, qubit + dim_space] stabiliser_matrix[:, [qubit, qubit + dim_space]] = stabiliser_matrix[:, [qubit + dim_space, qubit]] gates_list.append(gates.H(qubit)) qubits.append(qubit) @@ -358,12 +358,11 @@ def _col_reduce_x_matrix(stabiliser_matrix: np.ndarray, phases: np.ndarray) -> l phase_changes = stabiliser_matrix[:, col] + stabiliser_matrix[:, pivot_col + dim_space] + 1 phase_changes %= 2 phase_changes *= stabiliser_matrix[:, pivot_col] * stabiliser_matrix[:, col + dim_space] - phases += phase_changes + phases ^= phase_changes # X matrix: Add pivot column to column with 1 stabiliser_matrix[:, col] += stabiliser_matrix[:, pivot_col] # Z matrix: Add (column with 1)^th column to pivot column - stabiliser_matrix[:, pivot_col + dim_space] += stabiliser_matrix[:, col + dim_space] - stabiliser_matrix %= 2 + stabiliser_matrix[:, pivot_col + dim_space] ^= stabiliser_matrix[:, col + dim_space] gates_list.append(gates.CNOT(pivot_col, col)) pivot_col += 1 @@ -387,7 +386,7 @@ def _zero_z_matrix(stabiliser_matrix: np.ndarray, phases: np.ndarray) -> list[ga for i in range(dim): if stabiliser_matrix[i, dim_space + i] == 1: # For S(a)/H(a): r_i := r_i + x_{i,a} z_{i,a} for all i - phases += stabiliser_matrix[:, i] * stabiliser_matrix[:, i + dim_space] + phases ^= stabiliser_matrix[:, i] * stabiliser_matrix[:, i + dim_space] stabiliser_matrix[i, dim_space + i] = 0 s_gates.append(gates.S(i)) # Then remove the off-diagonal terms in each row From 680931f6bb3f0b0619cfaaf55395701dab2f3695 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 31 Aug 2026 14:04:10 +0800 Subject: [PATCH 33/36] Update other operations as well --- src/qibochem/measurement/util.py | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index 34384ff8..754c6f3a 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -120,7 +120,7 @@ def _binary_gaussian_elimination(vector_space: np.ndarray) -> np.ndarray: Performs Gaussian elimination on a binary vector_space. Returns the (unique) reduced row echelon form, and removes any zero rows as well """ - vector_space = np.array(vector_space) # Create a copy for returning + vector_space = np.array(vector_space, dtype=np.uint8) # Create a copy for returning rows, cols = vector_space.shape pivot_row = 0 @@ -141,8 +141,7 @@ def _binary_gaussian_elimination(vector_space: np.ndarray) -> np.ndarray: rows_to_reduce = rows_to_reduce[rows_to_reduce != pivot_row] # In GF(2), elimination is XOR with the pivot row. - vector_space[rows_to_reduce] += vector_space[pivot_row] - vector_space %= 2 + vector_space[rows_to_reduce] ^= vector_space[pivot_row] pivot_row += 1 if pivot_row == rows: @@ -158,7 +157,7 @@ def _binary_nullspace(binary_matrix: np.ndarray) -> np.ndarray: """Finds the nullspace of a binary_matrix, i.e. x s.t. Ax = 0""" dim = binary_matrix.shape[0] # Form the augmented matrix - aug_matrix = np.concatenate((binary_matrix.T, np.identity(binary_matrix.shape[1])), axis=1) + aug_matrix = np.concatenate((binary_matrix.T, np.identity(binary_matrix.shape[1], dtype=np.uint8)), axis=1) rref_aug_matrix = _binary_gaussian_elimination(aug_matrix) nullspace = rref_aug_matrix[dim:, dim:] return nullspace.astype(int) @@ -183,10 +182,9 @@ def _lagrangian_subspace(vector_space: np.ndarray) -> np.ndarray: space_to_orthogonalize = np.delete(vector_space, anticommuting_vector_indices, axis=0) for i1, vector in enumerate(space_to_orthogonalize): for i2, anticommuting_vector in enumerate(anticommuting_vectors): - space_to_orthogonalize[i1] += ( + space_to_orthogonalize[i1] ^= ( _symplectic_inner_product(vector, anticommuting_vectors[1 - i2]) * anticommuting_vector ) - space_to_orthogonalize = space_to_orthogonalize % 2 # Preferentially select Z over X first_nonzero_col = np.argmax(anticommuting_vectors, axis=1) @@ -233,18 +231,13 @@ def _get_sigma_terms(tau_terms: np.ndarray) -> tuple[np.ndarray, np.ndarray]: sigma_i = np.ravel(np.array([(0, 0) if _j != _i else _sigma_i for _j in range(dim)]).T) sigma_terms.append(sigma_i) # Orthogonalise the non-i^th terms: - new_tau_terms += np.array( + new_tau_terms ^= np.array( [ - # Not sure if need _j != _i or if _j > _i is good enough? - # Paper says do _j > _i, but then will have some non-commuting tau/sigma's...? _symplectic_inner_product(new_tau_terms[_j], sigma_i) * tau_i if _j != _i else np.zeros(2 * dim) - # symplectic_inner_product(new_tau_terms[_j], sigma_i) * tau_i if _j > _i else np.zeros(2 * dim) for _j in range(dim) ], dtype=np.uint8, ) - new_tau_terms = new_tau_terms % 2 - return new_tau_terms, np.array(sigma_terms, dtype=np.uint8) @@ -355,12 +348,14 @@ def _col_reduce_x_matrix(stabiliser_matrix: np.ndarray, phases: np.ndarray) -> l # Remove all nonzero entries on row _i using CNOT gates for col in nonzero_cols: # For CNOT(a, b): r_i := r_i + x_{i,a} z_{i,b} (x_{i,b} + z_{i,a} + 1), for all i - phase_changes = stabiliser_matrix[:, col] + stabiliser_matrix[:, pivot_col + dim_space] + 1 - phase_changes %= 2 - phase_changes *= stabiliser_matrix[:, pivot_col] * stabiliser_matrix[:, col + dim_space] + phase_changes = ( + stabiliser_matrix[:, pivot_col] + & stabiliser_matrix[:, col + dim_space] + & (stabiliser_matrix[:, col] ^ stabiliser_matrix[:, pivot_col + dim_space] ^ 1) + ) phases ^= phase_changes # X matrix: Add pivot column to column with 1 - stabiliser_matrix[:, col] += stabiliser_matrix[:, pivot_col] + stabiliser_matrix[:, col] ^= stabiliser_matrix[:, pivot_col] # Z matrix: Add (column with 1)^th column to pivot column stabiliser_matrix[:, pivot_col + dim_space] ^= stabiliser_matrix[:, col + dim_space] gates_list.append(gates.CNOT(pivot_col, col)) From 26d3a015437dd28986b60eabbe481f7f78b7efc0 Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 31 Aug 2026 14:54:41 +0800 Subject: [PATCH 34/36] Revert physical SWAPs after circuit synthesis --- src/qibochem/measurement/optimization.py | 8 ++------ tests/test_measurement_optimisation.py | 5 ----- 2 files changed, 2 insertions(+), 11 deletions(-) diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index 88ba0725..69b22b21 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -125,7 +125,6 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl ] # Otherwise, expression is a sum of terms term_list = [_term_to_string(term) for term in expression.args if _term_to_string(term)[0] in ("X", "Y", "Z")] - term_qubits = sorted({_get_qubit(op) for term in term_list for op in term.split()}) v_subspace = np.array([_pauli_to_symplectic(terms.split(), nqubits) for terms in term_list], dtype=np.uint8) v_basis = _binary_gaussian_elimination(v_subspace) @@ -143,25 +142,22 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl if method == "chong": x_result = _solve_linear_system(v_basis, v_subspace) # Map the solution onto the original set of qubits - qubit_map = dict(zip(sorted({q for pauli_op in x_result for q in pauli_op}), term_qubits)) phase_factors = [_phase_factor(v_basis[pauli_op]) for pauli_op in x_result] u_gates, phases = _synthesise_circuit(v_basis) - u_gates += [gates.SWAP(i, j) for i, j in qubit_map.items() if i != j] mapping = { - term: phase * prod(phases[i] * Z(qubit_map[i]) for i in soln) + term: phase * prod(phases[i] * Z(i) for i in soln) for term, phase, soln in zip(term_list, phase_factors, x_result) } elif method == "izmaylov": v_basis = _sort_tau_terms(v_basis) new_tau_terms, sigma_terms = _get_sigma_terms(v_basis) x_result = _solve_linear_system(new_tau_terms, v_subspace) - qubit_map = dict(zip(sorted({q for pauli_op in x_result for q in pauli_op}), term_qubits)) phase_factors = [_phase_factor(new_tau_terms[pauli_op]) for pauli_op in x_result] tau_term_str = [_symplectic_to_pauli(tau_i) for tau_i in new_tau_terms] sigma_term_str = [_symplectic_to_pauli(sigma_i) for sigma_i in sigma_terms] qwc_terms = [_symplectic_to_pauli(sum(sigma_terms[_x] for _x in pauli_op)) for pauli_op in x_result] mapping = { - term: phase * prod([getattr(symbols, sigma[0])(qubit_map[int(sigma[1:])]) for sigma in pauli_op]) + term: phase * prod([getattr(symbols, sigma[0])(int(sigma[1:])) for sigma in pauli_op]) for term, phase, pauli_op in zip(term_list, phase_factors, qwc_terms) } # Define the measurement gates diff --git a/tests/test_measurement_optimisation.py b/tests/test_measurement_optimisation.py index 4e21060f..b0dbbc77 100644 --- a/tests/test_measurement_optimisation.py +++ b/tests/test_measurement_optimisation.py @@ -18,11 +18,6 @@ def test_gc_measurement_mapping(): """Remaining coverage tests for _gc_measurement_mapping""" - # Check that qubits to measure remain unchanged - ham = SymbolicHamiltonian(X(2) * X(3) + Y(2) * Y(3) + Z(2) * Z(3), nqubits=4) - mapping, m_gates = _gc_measurement_mapping(ham.form, ham.nqubits, "chong") - assert {term.target_qubit for term in mapping.values() if hasattr(term, "target_qubit")} == {2, 3} - # Single term Hamiltonian ham = SymbolicHamiltonian(Z(2)) mapping, m_gates = _gc_measurement_mapping(ham.form, 2, "izmaylov") From 0496716e695711817454421d2208bb21e9a71eef Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 31 Aug 2026 14:56:50 +0800 Subject: [PATCH 35/36] Add test for earlier phase bug --- tests/test_expectation_samples.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_expectation_samples.py b/tests/test_expectation_samples.py index dc0f520a..00407295 100644 --- a/tests/test_expectation_samples.py +++ b/tests/test_expectation_samples.py @@ -108,6 +108,7 @@ def test_measurement_grouping_functionality(grouping, terms): ), (Y(1) * Y(2) + X(0) * X(1) * Z(2), 3, (gates.H(0), gates.H(1))), (Y(0) * X(1) + X(0) * Y(1) * Z(2), 3, (gates.H(0), gates.S(0), gates.H(1), gates.X(2), gates.H(2))), + (X(0) * X(1) * Z(2) * X(3) + X(0) * Y(2) * Y(3), 4, (gates.H(0), gates.H(1), gates.H(3))), ], ) def test_measurement_grouping_extra_tests(grouping, terms, nqubits, gates_to_add): From 15e384b77786110258bb821972c7739d7800bbbd Mon Sep 17 00:00:00 2001 From: Ah Wong <70616433+chmwzc@users.noreply.github.com> Date: Mon, 31 Aug 2026 15:36:50 +0800 Subject: [PATCH 36/36] Minor cleanup --- src/qibochem/measurement/optimization.py | 3 +-- src/qibochem/measurement/util.py | 2 +- tests/test_measurement_util.py | 1 - 3 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index 69b22b21..e745e6aa 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -16,7 +16,6 @@ from qibochem.measurement.util import ( _binary_gaussian_elimination, _binary_nullspace, - _get_qubit, _get_sigma_terms, _group_commuting_terms, _lagrangian_subspace, @@ -111,7 +110,7 @@ def _gc_measurement_mapping(expression: Expr, nqubits: int, method: str) -> tupl can be used to calculate the expectation values of ALL the terms in expression directly. Args: - expression (sympy.Expr): Group of Pauli terms that all mutually commute with each other qubitwise + expression (sympy.Expr): Group of Pauli terms that mutually commutes with each other nqubits (int): Number of qubits of the original Hamiltonian method (str): Circuit formulation to use, either "chong" (default) or "izmaylov" diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index 754c6f3a..1231a374 100644 --- a/src/qibochem/measurement/util.py +++ b/src/qibochem/measurement/util.py @@ -160,7 +160,7 @@ def _binary_nullspace(binary_matrix: np.ndarray) -> np.ndarray: aug_matrix = np.concatenate((binary_matrix.T, np.identity(binary_matrix.shape[1], dtype=np.uint8)), axis=1) rref_aug_matrix = _binary_gaussian_elimination(aug_matrix) nullspace = rref_aug_matrix[dim:, dim:] - return nullspace.astype(int) + return nullspace def _lagrangian_subspace(vector_space: np.ndarray) -> np.ndarray: diff --git a/tests/test_measurement_util.py b/tests/test_measurement_util.py index fe743a72..74fb48af 100644 --- a/tests/test_measurement_util.py +++ b/tests/test_measurement_util.py @@ -291,7 +291,6 @@ def test_zero_z_matrix(): ) phases = np.array([0, 0, 0, 0], dtype=np.uint8) gates_list = _zero_z_matrix(stabiliser_matrix, phases) - print(phases) # Single column operation, should have only CNOT gate assert len(gates_list) == 1 and gates_list[0].name == "s" assert phases[0] == 1