diff --git a/src/qibochem/measurement/optimization.py b/src/qibochem/measurement/optimization.py index cfa7acf5..e745e6aa 100644 --- a/src/qibochem/measurement/optimization.py +++ b/src/qibochem/measurement/optimization.py @@ -110,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" @@ -124,7 +124,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")] - 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] @@ -133,15 +133,19 @@ 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) + # 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 phase_factors = [_phase_factor(v_basis[pauli_op]) for pauli_op in x_result] - u_gates = _synthesise_circuit(v_basis) + u_gates, phases = _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(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) diff --git a/src/qibochem/measurement/util.py b/src/qibochem/measurement/util.py index b559c131..1231a374 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,108 +117,99 @@ 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 + 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) + vector_space = np.array(vector_space, dtype=np.uint8) # 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 - 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): - break + row = pivot_row + pivot_candidates[0] - # Always take the first nonzero column to sort - nonzero_cols = np.nonzero(np.any(subspace_to_sort, axis=0))[0] - _col = nonzero_cols[0] + # Swap current row with pivot row if needed. + if pivot_row != row: + vector_space[[row, pivot_row]] = vector_space[[pivot_row, row]] - col_indices = subspace_to_sort[:, _col].argsort()[::-1] - subspace_to_sort[:, :] = subspace_to_sort[col_indices] + # 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] - # 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] + # In GF(2), elimination is XOR with the pivot row. + vector_space[rows_to_reduce] ^= vector_space[pivot_row] - # 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 + pivot_row += 1 + if pivot_row == rows: + break # 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 + 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: """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) + return nullspace 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 - cp_vector_space = np.append([anticommuting_vectors[0]], space_to_orthogonalize, axis=0) - cp_vector_space = _binary_gaussian_elimination(cp_vector_space) + # Preferentially select Z over X + first_nonzero_col = np.argmax(anticommuting_vectors, axis=1) + selected_vector = anticommuting_vectors[np.argmax(first_nonzero_col)] + + 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: - """ - 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']] - """ - # 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 + """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 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 = { + 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]: @@ -228,7 +221,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] @@ -238,18 +231,14 @@ 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) - ] - ).astype(int) - new_tau_terms = new_tau_terms % 2 - - return new_tau_terms, np.array(sigma_terms) + ], + dtype=np.uint8, + ) + 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]: @@ -258,7 +247,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: @@ -293,9 +282,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: @@ -306,107 +295,126 @@ 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) + qubits = [] + zero_row_indices = np.where(np.all(x_matrix == 0, axis=1))[0] + 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 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) + break + zero_row_indices = np.where(np.all(x_matrix == 0, axis=1))[0] 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 """ 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: + # 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[:, 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] + # 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)) + 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. + 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 - 2. CZ gates used to remove off-diagonal entries on Z matrix + 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 """ 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? + 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 phases here + 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) + 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) + 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) - # 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 diff --git a/tests/test_expectation_samples.py b/tests/test_expectation_samples.py index 4155f949..00407295 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 @@ -31,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)) @@ -81,35 +66,58 @@ 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) + 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", + [ + (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), + 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))), + (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): + """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(): @@ -137,18 +145,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) @@ -156,24 +165,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, diff --git a/tests/test_measurement_optimisation.py b/tests/test_measurement_optimisation.py index c4f6a0cc..b0dbbc77 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 @@ -19,6 +18,7 @@ def test_gc_measurement_mapping(): """Remaining coverage tests for _gc_measurement_mapping""" + # 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 diff --git a/tests/test_measurement_util.py b/tests/test_measurement_util.py index 529cdcb8..74fb48af 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 @@ -168,15 +191,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 = [_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"]] + 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 @@ -231,11 +254,29 @@ 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) + 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 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( + [ + [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(): @@ -245,8 +286,11 @@ 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) + phases = np.array([0, 0, 0, 0], dtype=np.uint8) + gates_list = _zero_z_matrix(stabiliser_matrix, phases) # 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" + assert phases[0] == 1 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":