From d205fe2c31d0ad9875f45864fae4d64325087ce5 Mon Sep 17 00:00:00 2001 From: positr0nium Date: Thu, 7 Nov 2024 09:48:30 +0100 Subject: [PATCH 01/12] implemented prototypical ladder operator measurement --- src/qrisp/operators/qubit/measurement.py | 49 +++++++++++++-- .../operators/qubit/pauli_measurement.py | 60 ++++++++++++------- src/qrisp/operators/qubit/qubit_operator.py | 4 ++ src/qrisp/operators/qubit/qubit_term.py | 34 +++++++++++ src/qrisp/simulator/simulator.py | 2 - 5 files changed, 122 insertions(+), 27 deletions(-) diff --git a/src/qrisp/operators/qubit/measurement.py b/src/qrisp/operators/qubit/measurement.py index 0d45ccbf..8fd3aac2 100644 --- a/src/qrisp/operators/qubit/measurement.py +++ b/src/qrisp/operators/qubit/measurement.py @@ -140,7 +140,6 @@ def get_measurement( pauli_measurement = _measurement meas_circs = pauli_measurement.circuits - meas_qubits = pauli_measurement.qubits meas_ops = pauli_measurement.operators_int meas_coeffs = pauli_measurement.coefficients meas_shots = pauli_measurement.shots @@ -156,9 +155,10 @@ def get_measurement( for index,circuit in enumerate(meas_circs): curr = qc.copy() - curr.append(circuit.to_gate(), meas_qubits[index]) - - res = get_measurement_from_qc(curr, meas_qubits[index], backend, meas_shots[index]) + qubits = [qarg[i] for i in range(circuit.num_qubits())] + curr.append(circuit.to_gate(), qubits) + + res = get_measurement_from_qc(curr.transpile(), list(qarg), backend, meas_shots[index]) results.append(res) expectation = evaluate_expectation(results, meas_ops, meas_coeffs) @@ -199,6 +199,47 @@ def evaluate_observable(observable: int, x: int): else: return -1 +def evaluate_observable(observable: tuple, x: int): + """ + This method evaluates an observable that is a tensor product of Pauli-:math:`Z` operators + with respect to a measurement outcome. + + A Pauli operator of the form :math:`\prod_{i\in I}Z_i`, for some finite set of indices :math:`I\subset \mathbb N`, + is identified with an integer: + We identify the Pauli operator with the binary string that has ones at positions :math:`i\in I` + and zeros otherwise, and then convert this binary string to an integer. + + Parameters + ---------- + + observable : int + The observable represented as integer. + x : int + The measurement outcome represented as integer. + + Returns + ------- + int + The value of the observable with respect to the measurement outcome. + + """ + + z_int, AND_bits, AND_ctrl_state = observable + + sign_flip = bin(z_int & x).count('1') + + temp = (x ^ AND_ctrl_state) + + if AND_bits == 0: + return (-1)**sign_flip + + if temp & AND_bits == 0 or AND_bits == 0: + return (-1)**sign_flip/2 + else: + return 0 + + return (-1)**(sign_flip) + def evaluate_expectation(results, operators, coefficients): """ diff --git a/src/qrisp/operators/qubit/pauli_measurement.py b/src/qrisp/operators/qubit/pauli_measurement.py index 64e1a1c5..6e8a06ec 100644 --- a/src/qrisp/operators/qubit/pauli_measurement.py +++ b/src/qrisp/operators/qubit/pauli_measurement.py @@ -38,7 +38,7 @@ def __init__(self, operator): Parameters ---------- - operator: QubitOperator or BoundQubitOperator + operator: QubitHamiltonian or BoundQubitHamiltonian Attributes ---------- @@ -64,7 +64,7 @@ def __init__(self, operator): """ self.bases, self.operators_ind, self.operators_int, self.coefficients = self.commuting_qw_measurement(operator) self.variances, self.shots = self.measurement_shots() - self.circuits, self.qubits = self.measurement_circuits() + self.circuits = self.measurement_circuits() def commuting_qw_measurement(self, operator): """ @@ -72,6 +72,7 @@ def commuting_qw_measurement(self, operator): """ groups, bases = operator.commuting_qw_groups(show_bases=True) + self.groups = groups operators_ind = [] operators_int = [] coefficients = [] @@ -89,11 +90,11 @@ def commuting_qw_measurement(self, operator): curr_int = [] curr_coeff = [] - for factor,coeff in groups[i].terms_dict.items(): - ind = list(factor.factor_dict.keys()) + for term,coeff in groups[i].terms_dict.items(): + ind = list(term.factor_dict.keys()) curr_ind.append(ind) - curr_int.append(get_integer_from_indices(ind,positions[i])) + curr_int.append(term.serialize()) curr_coeff.append(float(coeff.real)) operators_ind.append(curr_ind) @@ -111,7 +112,8 @@ def measurement_shots(self): variances = [] for i in range(n): m = len(self.coefficients[i]) - var = sum(self.coefficients[i][j]**2 for j in range(m) if self.operators_int[i][j]>0) # Exclude constant term + + var = sum(self.coefficients[i][j]**2 for j in range(m) if (self.operators_int[i][j][0]+self.operators_int[i][j][1])>0) # Exclude constant term variances.append(var) N = sum(np.sqrt(x) for x in variances) shots = [np.sqrt(x)*N for x in variances] @@ -124,25 +126,41 @@ def measurement_circuits(self): """ circuits = [] - qubits = [] # Construct change of basis circuits - for basis in self.bases: - - basis_ = sorted(basis.factor_dict.items()) - qubits_ = sorted(basis.factor_dict.keys()) - - n = len(basis_) + for i in range(len(self.bases)): + + group = self.groups[i] + n = group.find_minimal_qubit_amount() qc = QuantumCircuit(n) - for i in range(n): - if basis_[i][1]=="X": - qc.h(i) - if basis_[i][1]=="Y": - qc.sx(i) + + factor_dict = self.bases[i].factor_dict + for j in range(n): + if factor_dict.get(j, "I")=="X": + qc.h(j) + if factor_dict.get(j, "I")=="Y": + qc.sx(j) + + for term in self.groups[i].terms_dict.keys(): + + ladder_operators = [base for base in term.factor_dict.items() if base[1] in ["A", "C"]] + + if len(ladder_operators): + anchor_factor = ladder_operators[-1] + + if anchor_factor[1] == "A": + qc.x(anchor_factor[0]) + + for j in range(len(ladder_operators)-1): + qc.cx(anchor_factor[0], j) + + if len(ladder_operators): + if anchor_factor[1] == "A": + qc.x(anchor_factor[0]) + + qc.h(anchor_factor[0]) circuits.append(qc) - qubits.append(qubits_) - return circuits, qubits + return circuits - diff --git a/src/qrisp/operators/qubit/qubit_operator.py b/src/qrisp/operators/qubit/qubit_operator.py index dbdeb08a..c7aae217 100644 --- a/src/qrisp/operators/qubit/qubit_operator.py +++ b/src/qrisp/operators/qubit/qubit_operator.py @@ -426,6 +426,10 @@ def subs(self, subs_dict): # Miscellaneous # + def find_minimal_qubit_amount(self): + indices = sum([list(term.factor_dict.keys()) for term in self.terms_dict.keys()], []) + return max(indices)+1 + def commutator(self, other): res = 0 diff --git a/src/qrisp/operators/qubit/qubit_term.py b/src/qrisp/operators/qubit/qubit_term.py index 33e57b70..5c1a8131 100644 --- a/src/qrisp/operators/qubit/qubit_term.py +++ b/src/qrisp/operators/qubit/qubit_term.py @@ -52,6 +52,40 @@ def copy(self): def is_identity(self): return len(self.factor_dict)==0 + def serialize(self): + z_int = 0 + and_int = 0 + ctrl_int = 0 + last_ladder_factor = None + factor_dict = self.factor_dict + for i in factor_dict.keys(): + bit = (1<= 1000000, no samples will be drawn and the distribution will #be returned instead @@ -200,7 +199,6 @@ def run(qc, shots, token="", iqs=None, insert_reset=True): #Generate samples else: - from numpy.random import choice # p_array = np.array(list(prob_dict.values())) From 8fe7a12fddd63b28482847279e4d9b04e7f9b4bc Mon Sep 17 00:00:00 2001 From: positr0nium Date: Thu, 7 Nov 2024 10:14:29 +0100 Subject: [PATCH 02/12] implemented to_pauli method --- .../reference/Hamiltonians/QubitOperator.rst | 10 +--- ....operators.qubit.QubitOperator.__add__.rst | 6 --- ...operators.qubit.QubitOperator.__iadd__.rst | 6 --- ...operators.qubit.QubitOperator.__imul__.rst | 6 --- ...operators.qubit.QubitOperator.__isub__.rst | 6 --- ....operators.qubit.QubitOperator.__mul__.rst | 6 --- ....operators.qubit.QubitOperator.__str__.rst | 6 --- ....operators.qubit.QubitOperator.__sub__.rst | 6 --- ...operators.qubit.QubitOperator.to_pauli.rst | 6 +++ src/qrisp/operators/qubit/qubit_operator.py | 48 +++++++++++++++++++ src/qrisp/operators/qubit/qubit_term.py | 23 +++++++++ .../test_to_pauli_converter.py | 34 +++++++++++++ 12 files changed, 113 insertions(+), 50 deletions(-) delete mode 100644 documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__add__.rst delete mode 100644 documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__iadd__.rst delete mode 100644 documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__imul__.rst delete mode 100644 documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__isub__.rst delete mode 100644 documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__mul__.rst delete mode 100644 documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__str__.rst delete mode 100644 documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__sub__.rst create mode 100644 documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.to_pauli.rst create mode 100644 tests/hamiltonian_tests/test_to_pauli_converter.py diff --git a/documentation/source/reference/Hamiltonians/QubitOperator.rst b/documentation/source/reference/Hamiltonians/QubitOperator.rst index 61ad3076..54b79363 100644 --- a/documentation/source/reference/Hamiltonians/QubitOperator.rst +++ b/documentation/source/reference/Hamiltonians/QubitOperator.rst @@ -17,15 +17,9 @@ Methods .. autosummary:: :toctree: generated/ - QubitOperator.__add__ - QubitOperator.__sub__ - QubitOperator.__mul__ - QubitOperator.__iadd__ - QubitOperator.__isub__ - QubitOperator.__imul__ - QubitOperator.__str__ QubitOperator.ground_state_energy QubitOperator.get_measurement QubitOperator.commuting_groups QubitOperator.commuting_qw_groups - QubitOperator.trotterization \ No newline at end of file + QubitOperator.trotterization + QubitOperator.to_pauli \ No newline at end of file diff --git a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__add__.rst b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__add__.rst deleted file mode 100644 index 4f40107f..00000000 --- a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__add__.rst +++ /dev/null @@ -1,6 +0,0 @@ -qrisp.operators.qubit.QubitOperator.\_\_add\_\_ -=============================================== - -.. currentmodule:: qrisp.operators.qubit - -.. automethod:: QubitOperator.__add__ \ No newline at end of file diff --git a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__iadd__.rst b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__iadd__.rst deleted file mode 100644 index 81b74aae..00000000 --- a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__iadd__.rst +++ /dev/null @@ -1,6 +0,0 @@ -qrisp.operators.qubit.QubitOperator.\_\_iadd\_\_ -================================================ - -.. currentmodule:: qrisp.operators.qubit - -.. automethod:: QubitOperator.__iadd__ \ No newline at end of file diff --git a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__imul__.rst b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__imul__.rst deleted file mode 100644 index a5502d9f..00000000 --- a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__imul__.rst +++ /dev/null @@ -1,6 +0,0 @@ -qrisp.operators.qubit.QubitOperator.\_\_imul\_\_ -================================================ - -.. currentmodule:: qrisp.operators.qubit - -.. automethod:: QubitOperator.__imul__ \ No newline at end of file diff --git a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__isub__.rst b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__isub__.rst deleted file mode 100644 index e719ee43..00000000 --- a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__isub__.rst +++ /dev/null @@ -1,6 +0,0 @@ -qrisp.operators.qubit.QubitOperator.\_\_isub\_\_ -================================================ - -.. currentmodule:: qrisp.operators.qubit - -.. automethod:: QubitOperator.__isub__ \ No newline at end of file diff --git a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__mul__.rst b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__mul__.rst deleted file mode 100644 index ba1cd58c..00000000 --- a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__mul__.rst +++ /dev/null @@ -1,6 +0,0 @@ -qrisp.operators.qubit.QubitOperator.\_\_mul\_\_ -=============================================== - -.. currentmodule:: qrisp.operators.qubit - -.. automethod:: QubitOperator.__mul__ \ No newline at end of file diff --git a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__str__.rst b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__str__.rst deleted file mode 100644 index 6edf1c99..00000000 --- a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__str__.rst +++ /dev/null @@ -1,6 +0,0 @@ -qrisp.operators.qubit.QubitOperator.\_\_str\_\_ -=============================================== - -.. currentmodule:: qrisp.operators.qubit - -.. automethod:: QubitOperator.__str__ \ No newline at end of file diff --git a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__sub__.rst b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__sub__.rst deleted file mode 100644 index 4a4b9d37..00000000 --- a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.__sub__.rst +++ /dev/null @@ -1,6 +0,0 @@ -qrisp.operators.qubit.QubitOperator.\_\_sub\_\_ -=============================================== - -.. currentmodule:: qrisp.operators.qubit - -.. automethod:: QubitOperator.__sub__ \ No newline at end of file diff --git a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.to_pauli.rst b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.to_pauli.rst new file mode 100644 index 00000000..3242b676 --- /dev/null +++ b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.to_pauli.rst @@ -0,0 +1,6 @@ +qrisp.operators.qubit.QubitOperator.to\_pauli +============================================= + +.. currentmodule:: qrisp.operators.qubit + +.. automethod:: QubitOperator.to_pauli \ No newline at end of file diff --git a/src/qrisp/operators/qubit/qubit_operator.py b/src/qrisp/operators/qubit/qubit_operator.py index c7aae217..376e8534 100644 --- a/src/qrisp/operators/qubit/qubit_operator.py +++ b/src/qrisp/operators/qubit/qubit_operator.py @@ -537,6 +537,54 @@ def recursive_TP(keys,term_dict): # res.sum_duplicates() return M + def to_array(self, factor_amount = None): + """ + Returns a numpy array describing the operator. + + Parameters + ---------- + factor_amount : int, optional + The amount of factors to represent this array. The array will have + the dimension $2^n \times 2^n$, where n is the amount of factors. + By default the minimal number is chosen. + + Returns + ------- + np.ndarray + The array describing the operator. + + """ + return self.to_sparse_matrix(factor_amount).todense() + + def to_pauli(self): + """ + Returns an equivalent operator, which however only contains Pauli factors. + + Returns + ------- + QubitOperator + An operator that contains only Pauli-Factor. + + Examples + -------- + + We create a QubitOperator containing A and C terms and convert it to a + Pauli based representation. + + >>> from qrisp.operators import A,C,Z + >>> H = A(0)*C(1)*Z(2) + >>> print(H.to_pauli()) + 0.25*X_0*X_1*Z_2 + 0.25*I*X_0*Y_1*Z_2 - 0.25*I*Y_0*X_1*Z_2 + 0.25*Y_0*Y_1*Z_2 + + """ + + res = 0 + for term, coeff in self.terms_dict.items(): + res += coeff*term.to_pauli() + if isinstance(res, (float, int)): + return QubitOperator({QubitTerm({}) : res}) + return res + def ground_state_energy(self): """ Calculates the ground state energy (i.e., the minimum eigenvalue) of the operator classically. diff --git a/src/qrisp/operators/qubit/qubit_term.py b/src/qrisp/operators/qubit/qubit_term.py index 5c1a8131..164d8e26 100644 --- a/src/qrisp/operators/qubit/qubit_term.py +++ b/src/qrisp/operators/qubit/qubit_term.py @@ -86,6 +86,29 @@ def serialize(self): return (z_int, and_int, ctrl_int) + def to_pauli(self): + + from qrisp.operators import X, Y, Z + res = 1 + + for i, factor in self.factor_dict.items(): + if factor == "X": + res *= X(i) + elif factor == "Y": + res *= Y(i) + elif factor == "Z": + res *= Z(i) + elif factor == "A": + res *= (X(i) - 1j*Y(i))*0.5 + elif factor == "C": + res *= (X(i) + 1j*Y(i))*0.5 + elif factor == "P0": + res *= (Z(i) + 1)*0.5 + elif factor == "P1": + res *= (Z(i) - 1)*(-0.5) + + return res + # # Simulation # diff --git a/tests/hamiltonian_tests/test_to_pauli_converter.py b/tests/hamiltonian_tests/test_to_pauli_converter.py new file mode 100644 index 00000000..2badf045 --- /dev/null +++ b/tests/hamiltonian_tests/test_to_pauli_converter.py @@ -0,0 +1,34 @@ +""" +\******************************************************************************** +* Copyright (c) 2024 the Qrisp authors +* +* This program and the accompanying materials are made available under the +* terms of the Eclipse Public License 2.0 which is available at +* http://www.eclipse.org/legal/epl-2.0. +* +* This Source Code may also be made available under the following Secondary +* Licenses when the conditions for such availability set forth in the Eclipse +* Public License, v. 2.0 are satisfied: GNU General Public License, version 2 +* with the GNU Classpath Exception which is +* available at https://www.gnu.org/software/classpath/license.html. +* +* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0 +********************************************************************************/ +""" + +from qrisp.operators import X, Y, Z, A, C, P0, P1 +from numpy.linalg import norm + +def test_to_pauli_converter(): + + operator_list = [lambda x : 1, X, Y, Z, A, C, P0, P1] + + for O0 in operator_list: + for O1 in operator_list: + for O2 in operator_list: + H = O0(0)*O1(1)*O2(2) + if isinstance(H, int): + continue + + assert norm(H.to_array() - H.to_pauli().to_array()) < 1E-5 + From 6fbc1fc373d8045706e86d12381dbbf3be0b9178 Mon Sep 17 00:00:00 2001 From: positr0nium Date: Thu, 7 Nov 2024 11:49:05 +0100 Subject: [PATCH 03/12] wrote an exhaustive test for ladder operator measurement and fixed all appearing bugs --- .../logic_synthesis/truth_tables.py | 46 +++++++-------- src/qrisp/operators/qubit/measurement.py | 50 ++++++---------- .../operators/qubit/pauli_measurement.py | 2 +- src/qrisp/operators/qubit/qubit_term.py | 8 ++- .../test_measurement_method.py | 59 +++++++++++++++++++ 5 files changed, 105 insertions(+), 60 deletions(-) create mode 100644 tests/hamiltonian_tests/test_measurement_method.py diff --git a/src/qrisp/alg_primitives/logic_synthesis/truth_tables.py b/src/qrisp/alg_primitives/logic_synthesis/truth_tables.py index 29a2232f..4766f55b 100644 --- a/src/qrisp/alg_primitives/logic_synthesis/truth_tables.py +++ b/src/qrisp/alg_primitives/logic_synthesis/truth_tables.py @@ -339,25 +339,25 @@ def rw_spectrum(f): return np.dot(T(n), f) -def C(f): - size = len(f) +# def C(f): +# size = len(f) - if np.log2(size) != int(np.log2(size)): - raise Exception( - "The given function does not have the length to properly represent " - "a truth table" - ) +# if np.log2(size) != int(np.log2(size)): +# raise Exception( +# "The given function does not have the length to properly represent " +# "a truth table" +# ) - n = int(np.log2(size)) +# n = int(np.log2(size)) - rw_spec = rw_spectrum(f) - sum_ = 0 - for i in range(size): - sum_ += sum(int_as_array(i, n)) * rw_spec[i] ** 2 +# rw_spec = rw_spectrum(f) +# sum_ = 0 +# for i in range(size): +# sum_ += sum(int_as_array(i, n)) * rw_spec[i] ** 2 - sum_ = sum_ / 2 ** (n - 2) +# sum_ = sum_ / 2 ** (n - 2) - return 1 / 2 * (n * size - sum_) +# return 1 / 2 * (n * size - sum_) def NZ(f): @@ -370,17 +370,17 @@ def NZ(f): return sum_ -def D(f): - size = len(f) - if np.log2(size) != int(np.log2(size)): - raise Exception( - "The given function does not have the length to properly represent " - "a truth table" - ) +# def D(f): +# size = len(f) +# if np.log2(size) != int(np.log2(size)): +# raise Exception( +# "The given function does not have the length to properly represent " +# "a truth table" +# ) - n = int(np.log2(size)) +# n = int(np.log2(size)) - return int(n * 2 ** (n - 3) * NZ(f) + C(f)) +# return int(n * 2 ** (n - 3) * NZ(f) + C(f)) def synth_poly(truth_table, column=0, coeff=None): diff --git a/src/qrisp/operators/qubit/measurement.py b/src/qrisp/operators/qubit/measurement.py index 8fd3aac2..df2e34b3 100644 --- a/src/qrisp/operators/qubit/measurement.py +++ b/src/qrisp/operators/qubit/measurement.py @@ -158,6 +158,7 @@ def get_measurement( qubits = [qarg[i] for i in range(circuit.num_qubits())] curr.append(circuit.to_gate(), qubits) + # print(curr.transpile()) res = get_measurement_from_qc(curr.transpile(), list(qarg), backend, meas_shots[index]) results.append(res) @@ -169,35 +170,6 @@ def get_measurement( # Evaluate expectation # -def evaluate_observable(observable: int, x: int): - """ - This method evaluates an observable that is a tensor product of Pauli-:math:`Z` operators - with respect to a measurement outcome. - - A Pauli operator of the form :math:`\prod_{i\in I}Z_i`, for some finite set of indices :math:`I\subset \mathbb N`, - is identified with an integer: - We identify the Pauli operator with the binary string that has ones at positions :math:`i\in I` - and zeros otherwise, and then convert this binary string to an integer. - - Parameters - ---------- - - observable : int - The observable represented as integer. - x : int - The measurement outcome represented as integer. - - Returns - ------- - int - The value of the observable with respect to the measurement outcome. - - """ - - if bin(observable & x).count('1') % 2 == 0: - return 1 - else: - return -1 def evaluate_observable(observable: tuple, x: int): """ @@ -224,20 +196,32 @@ def evaluate_observable(observable: tuple, x: int): """ - z_int, AND_bits, AND_ctrl_state = observable + z_int, AND_bits, AND_ctrl_state, contains_ladder = observable sign_flip = bin(z_int & x).count('1') + from qrisp import bin_rep temp = (x ^ AND_ctrl_state) + # if AND_bits != 0: + # print(bin_rep(AND_ctrl_state, 3)) + # print(bin_rep(AND_bits, 3)) + # print(bin_rep(temp, 3)) + # print(bin_rep(z_int, 3)) + # print("=====") + if contains_ladder: + prefactor = 0.5 + else: + prefactor = 1 if AND_bits == 0: - return (-1)**sign_flip + return prefactor*(-1)**sign_flip - if temp & AND_bits == 0 or AND_bits == 0: - return (-1)**sign_flip/2 + if temp & AND_bits == 0: + return prefactor*(-1)**sign_flip else: return 0 + return (-1)**(sign_flip) diff --git a/src/qrisp/operators/qubit/pauli_measurement.py b/src/qrisp/operators/qubit/pauli_measurement.py index 6e8a06ec..797f78b9 100644 --- a/src/qrisp/operators/qubit/pauli_measurement.py +++ b/src/qrisp/operators/qubit/pauli_measurement.py @@ -152,7 +152,7 @@ def measurement_circuits(self): qc.x(anchor_factor[0]) for j in range(len(ladder_operators)-1): - qc.cx(anchor_factor[0], j) + qc.cx(anchor_factor[0], ladder_operators[j][0]) if len(ladder_operators): if anchor_factor[1] == "A": diff --git a/src/qrisp/operators/qubit/qubit_term.py b/src/qrisp/operators/qubit/qubit_term.py index 164d8e26..801c7e8b 100644 --- a/src/qrisp/operators/qubit/qubit_term.py +++ b/src/qrisp/operators/qubit/qubit_term.py @@ -63,13 +63,13 @@ def serialize(self): if factor_dict[i] in ["X", "Y", "Z"]: z_int |= bit - last_ladder_factor continue elif factor_dict[i] == "A": ctrl_int |= bit - last_ladder_factor + last_ladder_factor = bit pass elif factor_dict[i] == "C": + last_ladder_factor = bit pass elif factor_dict[i] == "P0": pass @@ -82,9 +82,11 @@ def serialize(self): and_int |= bit if last_ladder_factor is not None: + pass and_int ^= last_ladder_factor + z_int ^= last_ladder_factor - return (z_int, and_int, ctrl_int) + return (z_int, and_int, ctrl_int, last_ladder_factor is not None) def to_pauli(self): diff --git a/tests/hamiltonian_tests/test_measurement_method.py b/tests/hamiltonian_tests/test_measurement_method.py new file mode 100644 index 00000000..a00c8949 --- /dev/null +++ b/tests/hamiltonian_tests/test_measurement_method.py @@ -0,0 +1,59 @@ +""" +\******************************************************************************** +* Copyright (c) 2024 the Qrisp authors +* +* This program and the accompanying materials are made available under the +* terms of the Eclipse Public License 2.0 which is available at +* http://www.eclipse.org/legal/epl-2.0. +* +* This Source Code may also be made available under the following Secondary +* Licenses when the conditions for such availability set forth in the Eclipse +* Public License, v. 2.0 are satisfied: GNU General Public License, version 2 +* with the GNU Classpath Exception which is +* available at https://www.gnu.org/software/classpath/license.html. +* +* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0 +********************************************************************************/ +""" + +from qrisp.operators import X, Y, Z, A, C, P0, P1 +from numpy.linalg import norm +from qrisp import * + +def test_measurement_method(): + + def testing_helper(qv): + operator_list = [lambda x : 1, X, Y, Z, A, C, P0, P1] + for O0 in operator_list: + for O1 in operator_list: + for O2 in operator_list: + for O3 in operator_list: + H = O0(0)*O1(1)*O2(2)*O3(3) + if isinstance(H, int): + continue + + print(H) + assert abs(H.get_measurement(qv, precision = 0.001) - H.to_pauli().get_measurement(qv, precision = 0.001)) < 1E-2 + + qv = QuantumVariable(4) + + testing_helper(qv) + + h(qv[0]) + + testing_helper(qv) + + cx(qv[0], qv[1]) + + testing_helper(qv) + + cx(qv[0], qv[2]) + + testing_helper(qv) + + h(qv[0]) + + + + + \ No newline at end of file From 81cdb13011a1b2c03550bbf16d5f203f9193f111 Mon Sep 17 00:00:00 2001 From: positr0nium Date: Thu, 7 Nov 2024 12:58:52 +0100 Subject: [PATCH 04/12] implemented the QubitOperator.hermitize --- .../reference/Hamiltonians/QubitOperator.rst | 4 +- ....operators.qubit.QubitOperator.adjoint.rst | 6 +++ ...perators.qubit.QubitOperator.hermitize.rst | 6 +++ src/qrisp/operators/qubit/qubit_operator.py | 46 ++++++++++++++++++- src/qrisp/operators/qubit/qubit_term.py | 36 ++++++++++++--- 5 files changed, 89 insertions(+), 9 deletions(-) create mode 100644 documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.adjoint.rst create mode 100644 documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.hermitize.rst diff --git a/documentation/source/reference/Hamiltonians/QubitOperator.rst b/documentation/source/reference/Hamiltonians/QubitOperator.rst index 54b79363..48a64e3c 100644 --- a/documentation/source/reference/Hamiltonians/QubitOperator.rst +++ b/documentation/source/reference/Hamiltonians/QubitOperator.rst @@ -22,4 +22,6 @@ Methods QubitOperator.commuting_groups QubitOperator.commuting_qw_groups QubitOperator.trotterization - QubitOperator.to_pauli \ No newline at end of file + QubitOperator.to_pauli + QubitOperator.adjoint + QubitOperator.hermitize \ No newline at end of file diff --git a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.adjoint.rst b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.adjoint.rst new file mode 100644 index 00000000..53df69e1 --- /dev/null +++ b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.adjoint.rst @@ -0,0 +1,6 @@ +qrisp.operators.qubit.QubitOperator.adjoint +=========================================== + +.. currentmodule:: qrisp.operators.qubit + +.. automethod:: QubitOperator.adjoint \ No newline at end of file diff --git a/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.hermitize.rst b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.hermitize.rst new file mode 100644 index 00000000..ac2aea66 --- /dev/null +++ b/documentation/source/reference/Hamiltonians/generated/qrisp.operators.qubit.QubitOperator.hermitize.rst @@ -0,0 +1,6 @@ +qrisp.operators.qubit.QubitOperator.hermitize +============================================= + +.. currentmodule:: qrisp.operators.qubit + +.. automethod:: QubitOperator.hermitize \ No newline at end of file diff --git a/src/qrisp/operators/qubit/qubit_operator.py b/src/qrisp/operators/qubit/qubit_operator.py index 376e8534..f1082182 100644 --- a/src/qrisp/operators/qubit/qubit_operator.py +++ b/src/qrisp/operators/qubit/qubit_operator.py @@ -90,7 +90,7 @@ class QubitOperator(Hamiltonian): H = 1+2*X(0)+3*X(0)*Y(1)*A(2)+C(4)*P1(0) H - Yields $3*A_2*X_0*Y_1 + C_4*P1_0 + 1 + 2*X_0$. + Yields $1 + P^1_0*C_4 + 2*X_0 + 3*X_0*Y_1*A_2$. Investigate the simulation circuit by simulating for a symbolic amount of time: @@ -353,7 +353,7 @@ def __isub__(self,other): """ if isinstance(other,(int,float,complex)): - self.terms_dict[termTerm()] = self.terms_dict.get(termTerm(),0)-other + self.terms_dict[QubitTerm()] = self.terms_dict.get(QubitTerm(),0)-other return self if not isinstance(other,QubitOperator): raise TypeError("Cannot add QubitOperator and "+str(type(other))) @@ -585,6 +585,47 @@ def to_pauli(self): return QubitOperator({QubitTerm({}) : res}) return res + def adjoint(self): + """ + Returns an the adjoint operator. + + Returns + ------- + QubitOperator + The adjoint operator. + + Examples + -------- + + We create a QubitOperator and inspect it' adjoint. + + >>> from qrisp.operators import A,C,Z + >>> H = A(0)*C(1)*Z(2) + >>> print(H.adjoint()) + C_0*A_1*Z_2 + """ + res = 0 + for term, coeff in self.terms_dict.items(): + res += coeff*term.adjoint() + if isinstance(res, (float, int)): + return QubitOperator({QubitTerm({}) : res}) + return res + + def hermitize(self): + """ + Returns the hermitian part of self. + + $H = (O + O^\dagger)/2$ + + Returns + ------- + QubitOperator + The hermitian part. + + """ + return 0.5*(self + self.adjoint()) + + def ground_state_energy(self): """ Calculates the ground state energy (i.e., the minimum eigenvalue) of the operator classically. @@ -872,3 +913,4 @@ def U(qarg, t=1, steps=1, iter=1): trotter_step(qarg, t, steps) return U + diff --git a/src/qrisp/operators/qubit/qubit_term.py b/src/qrisp/operators/qubit/qubit_term.py index 801c7e8b..48d8c954 100644 --- a/src/qrisp/operators/qubit/qubit_term.py +++ b/src/qrisp/operators/qubit/qubit_term.py @@ -111,6 +111,28 @@ def to_pauli(self): return res + def adjoint(self): + from qrisp.operators import X, Y, Z, A, C, P0, P1 + res = 1 + + for i, factor in self.factor_dict.items(): + if factor == "X": + res *= X(i) + elif factor == "Y": + res *= Y(i) + elif factor == "Z": + res *= Z(i) + elif factor == "A": + res *= C(i) + elif factor == "C": + res *= A(i) + elif factor == "P0": + res *= P0(i) + elif factor == "P1": + res *= P1(i) + + return res + # # Simulation # @@ -397,17 +419,19 @@ def to_spin(P, index): if P=="Z": return Z_(index) if P=="A": - return Symbol("A_" + str(index)) + return Symbol("A_" + str(index), commutative = False) if P=="C": - return Symbol("C_" + str(index)) + return Symbol("C_" + str(index), commutative = False) if P=="P0": - return Symbol("P0_" + str(index)) + return Symbol("P^0_" + str(index), commutative = False) if P=="P1": - return Symbol("P1_" + str(index)) + return Symbol("P^1_" + str(index), commutative = False) expr = 1 - for index,P in self.factor_dict.items(): - expr *= to_spin(P,str(index)) + index_list = sorted(list(self.factor_dict.keys())) + + for i in index_list: + expr = expr*to_spin(self.factor_dict[i],str(i)) return expr From ac860b56da2f149bb44ee7efce690c73adeee335 Mon Sep 17 00:00:00 2001 From: positr0nium Date: Thu, 7 Nov 2024 16:08:08 +0100 Subject: [PATCH 05/12] made the new measurement procedure available for BoundQubitOperator --- .../operators/qubit/bound_qubit_operator.py | 46 +++++++++++++++---- src/qrisp/operators/qubit/bound_qubit_term.py | 4 +- src/qrisp/operators/qubit/measurement.py | 6 +-- src/qrisp/operators/qubit/qubit_operator.py | 2 +- tests/test_hamiltonian.py | 13 +++--- 5 files changed, 48 insertions(+), 23 deletions(-) diff --git a/src/qrisp/operators/qubit/bound_qubit_operator.py b/src/qrisp/operators/qubit/bound_qubit_operator.py index 38d4e967..fe3d37e8 100644 --- a/src/qrisp/operators/qubit/bound_qubit_operator.py +++ b/src/qrisp/operators/qubit/bound_qubit_operator.py @@ -644,20 +644,46 @@ def get_measurement( #Yields 1.0 """ - return get_measurement(self, - qarg, - precision=precision, - backend=backend, - shots=shots, - compile=compile, - compilation_kwargs=compilation_kwargs, - subs_dic=subs_dic, - precompiled_qc=precompiled_qc, - _measurement=_measurement) + + unbound_operator, qarg = self.unbind() + + return unbound_operator.get_measurement(qarg, + precision=precision, + backend=backend, + shots=shots, + compile=compile, + compilation_kwargs=compilation_kwargs, + subs_dic=subs_dic, + precompiled_qc=precompiled_qc, + _measurement=_measurement) + + def unbind(self): + from qrisp.operators import QubitTerm, QubitOperator + + participating_qubits = [] + for term, coeff in self.terms_dict.items(): + for qb in term.factor_dict.keys(): + if qb not in participating_qubits: + participating_qubits.append(qb) + + index_inv = {participating_qubits[i] : i for i in range(len(participating_qubits))} + + unbound_operator = 0 + + for term, coeff in self.terms_dict.items(): + new_term_factor_dict = {} + for qb, factor in term.factor_dict.items(): + new_term_factor_dict[index_inv[qb]] = factor + + unbound_operator += QubitOperator({QubitTerm(new_term_factor_dict): coeff}) + + return unbound_operator, participating_qubits + # # Trotterization # + def trotterization(self): r""" diff --git a/src/qrisp/operators/qubit/bound_qubit_term.py b/src/qrisp/operators/qubit/bound_qubit_term.py index 10b6d554..c86a93c9 100644 --- a/src/qrisp/operators/qubit/bound_qubit_term.py +++ b/src/qrisp/operators/qubit/bound_qubit_term.py @@ -205,6 +205,4 @@ def commute_qw(self, other): for key in keys: if a.get(key,"I")!="I" and b.get(key,"I")!="I" and a.get(key,"I")!=b.get(key,"I"): return False - return True - - + return True \ No newline at end of file diff --git a/src/qrisp/operators/qubit/measurement.py b/src/qrisp/operators/qubit/measurement.py index df2e34b3..bf315f20 100644 --- a/src/qrisp/operators/qubit/measurement.py +++ b/src/qrisp/operators/qubit/measurement.py @@ -94,12 +94,12 @@ def get_measurement( qs = qarg.qs elif isinstance(qarg,list): qs = QuantumSession() - for qv in qarg: - if qv.is_deleted(): + for arg in qarg: + if isinstance(arg, QuantumVariable) and qv.is_deleted(): raise Exception( "Tried to measure QuantumArray containing deleted QuantumVariables" ) - merge(qs,qv.qs) + merge(qs,arg) if backend is None: if qs.backend is None: diff --git a/src/qrisp/operators/qubit/qubit_operator.py b/src/qrisp/operators/qubit/qubit_operator.py index f1082182..adc92dd7 100644 --- a/src/qrisp/operators/qubit/qubit_operator.py +++ b/src/qrisp/operators/qubit/qubit_operator.py @@ -639,7 +639,7 @@ def ground_state_energy(self): from scipy.sparse.linalg import eigsh - M = self.to_sparse_matrix() + M = self.hermitize().to_sparse_matrix() # Compute the smallest eigenvalue eigenvalues, _ = eigsh(M, k=1, which='SA') # 'SA' stands for smallest algebraic E = eigenvalues[0] diff --git a/tests/test_hamiltonian.py b/tests/test_hamiltonian.py index 2175234b..2040f20d 100644 --- a/tests/test_hamiltonian.py +++ b/tests/test_hamiltonian.py @@ -28,12 +28,13 @@ def test_pauli_hamiltonian(): res = H.get_measurement(qv) assert np.abs(res-0.0) < 2e-2 - qtype = QuantumVariable(2) - q_array = QuantumArray(qtype, shape=(2)) - h(q_array) - H = Z(0)*Z(1) + X(2)*X(3) - res = H.get_measurement(q_array) - assert np.abs(res-1.0) < 2e-2 + # What is the semantics here? + # qtype = QuantumVariable(2) + # q_array = QuantumArray(qtype, shape=(2)) + # h(q_array) + # H = Z(0)*Z(1) + X(2)*X(3) + # res = H.get_measurement(q_array) + # assert np.abs(res-1.0) < 2e-2 def test_bound_pauli_hamiltonian(): From dd30ab2966a4d62c82159e0694ece5ead1561711 Mon Sep 17 00:00:00 2001 From: positr0nium Date: Fri, 8 Nov 2024 12:49:29 +0100 Subject: [PATCH 06/12] deprecated PauliMeasurement class --- src/qrisp/algorithms/vqe/vqe_problem.py | 23 ++--- src/qrisp/operators/qubit/__init__.py | 2 +- .../operators/qubit/bound_qubit_operator.py | 1 - src/qrisp/operators/qubit/measurement.py | 69 +++++++++------ ...i_measurement.py => pauli_measurement_.py} | 0 src/qrisp/operators/qubit/qubit_operator.py | 86 +++++++++++++++++-- src/qrisp/operators/qubit/qubit_term.py | 25 ++---- tests/{ => operator_tests}/test_VQE.py | 0 .../test_VQE_electronic_structure.py | 0 .../test_VQE_heisenberg.py | 0 .../test_fermionic_hamiltonian_simulation.py | 0 .../test_fermionic_term.py | 0 .../test_fermionic_to_qubit.py | 0 .../{ => operator_tests}/test_hamiltonian.py | 0 .../test_measurement_method.py | 2 +- .../test_qubit_hamiltonian_commutator.py | 0 .../test_qubit_hamiltonian_simulation.py | 0 .../test_to_pauli_converter.py | 0 18 files changed, 136 insertions(+), 72 deletions(-) rename src/qrisp/operators/qubit/{pauli_measurement.py => pauli_measurement_.py} (100%) rename tests/{ => operator_tests}/test_VQE.py (100%) rename tests/{ => operator_tests}/test_VQE_electronic_structure.py (100%) rename tests/{ => operator_tests}/test_VQE_heisenberg.py (100%) rename tests/{hamiltonian_tests => operator_tests}/test_fermionic_hamiltonian_simulation.py (100%) rename tests/{hamiltonian_tests => operator_tests}/test_fermionic_term.py (100%) rename tests/{hamiltonian_tests => operator_tests}/test_fermionic_to_qubit.py (100%) rename tests/{ => operator_tests}/test_hamiltonian.py (100%) rename tests/{hamiltonian_tests => operator_tests}/test_measurement_method.py (93%) rename tests/{hamiltonian_tests => operator_tests}/test_qubit_hamiltonian_commutator.py (100%) rename tests/{hamiltonian_tests => operator_tests}/test_qubit_hamiltonian_simulation.py (100%) rename tests/{hamiltonian_tests => operator_tests}/test_to_pauli_converter.py (100%) diff --git a/src/qrisp/algorithms/vqe/vqe_problem.py b/src/qrisp/algorithms/vqe/vqe_problem.py index 518483fa..17c76e5c 100644 --- a/src/qrisp/algorithms/vqe/vqe_problem.py +++ b/src/qrisp/algorithms/vqe/vqe_problem.py @@ -172,7 +172,7 @@ def compile_circuit(self, qarg, depth): return compiled_qc, theta - def optimization_routine(self, qarg, depth, mes_kwargs, measurement, max_iter, init_type="random", init_point=None, optimizer="COBYLA"): + def optimization_routine(self, qarg, depth, mes_kwargs, max_iter, init_type="random", init_point=None, optimizer="COBYLA"): """ Wrapper subroutine for the optimization method used in QAOA. The initial values are set and the optimization via ``COBYLA`` is conducted here. @@ -184,8 +184,6 @@ def optimization_routine(self, qarg, depth, mes_kwargs, measurement, max_iter, i The amont of VQE layers. mes_kwargs : dict The keyword arguments for the measurement function. - measurement : PauliMeasurement - The measurement setttings for the measurement function. max_iter : int The maximum number of iterations for the optimization method. init_type : string, optional @@ -204,7 +202,7 @@ def optimization_routine(self, qarg, depth, mes_kwargs, measurement, max_iter, i """ # Define optimization wrapper function to be minimized using VQE - def optimization_wrapper(theta, qc, symbols, qarg, mes_kwargs, measurement): + def optimization_wrapper(theta, qc, symbols, qarg, mes_kwargs): """ Wrapper function for the optimization method used in VQE. @@ -222,8 +220,6 @@ def optimization_wrapper(theta, qc, symbols, qarg, mes_kwargs, measurement): The duplicated quantum variable to which the quantum circuit is applied. mes_kwargs : dict The keyword arguments for the measurement function. - measurement : PauliMeasurement - The measurement setttings for the measurement function. Returns ------- @@ -233,7 +229,7 @@ def optimization_wrapper(theta, qc, symbols, qarg, mes_kwargs, measurement): subs_dic = {symbols[i] : theta[i] for i in range(len(symbols))} - expectation = self.hamiltonian.get_measurement(qarg, subs_dic = subs_dic, precompiled_qc = qc, _measurement=measurement, **mes_kwargs) + expectation = self.hamiltonian.get_measurement(qarg, subs_dic = subs_dic, precompiled_qc = qc, **mes_kwargs) if self.callback: self.optimization_costs.append(expectation) @@ -257,7 +253,7 @@ def optimization_wrapper(theta, qc, symbols, qarg, mes_kwargs, measurement): init_point, method=optimizer, options={'maxiter':max_iter}, - args = (compiled_qc, symbols, qarg, mes_kwargs, measurement)) + args = (compiled_qc, symbols, qarg, mes_kwargs)) return res_sample['x'] @@ -301,10 +297,8 @@ def run(self, qarg, depth, mes_kwargs = {}, max_iter = 50, init_type = "random", if not "precision" in mes_kwargs: mes_kwargs["precision"] = 0.01 - # Measurement settings - measurement = self.hamiltonian.pauli_measurement() - optimal_theta = self.optimization_routine(qarg, depth, mes_kwargs, measurement, max_iter, init_type, init_point, optimizer) + optimal_theta = self.optimization_routine(qarg, depth, mes_kwargs, max_iter, init_type, init_point, optimizer) # Prepare the initial state for particular problem instance, the default is the \ket{0} state if self.init_function is not None: @@ -314,7 +308,7 @@ def run(self, qarg, depth, mes_kwargs = {}, max_iter = 50, init_type = "random", for i in range(depth): self.ansatz_function(qarg,[optimal_theta[self.num_params*i+j] for j in range(self.num_params)]) - opt_res = self.hamiltonian.get_measurement(qarg,_measurement=measurement,**mes_kwargs) + opt_res = self.hamiltonian.get_measurement(qarg,**mes_kwargs) return opt_res @@ -356,10 +350,7 @@ def train_function(self, qarg, depth, mes_kwargs = {}, max_iter = 50, init_type if not "precision" in mes_kwargs: mes_kwargs["precision"] = 0.01 - # Measurement settings - measurement = self.hamiltonian.pauli_measurement() - - optimal_theta = self.optimization_routine(qarg, depth, mes_kwargs, measurement, max_iter, init_type, init_point, optimizer) + optimal_theta = self.optimization_routine(qarg, depth, mes_kwargs, max_iter, init_type, init_point, optimizer) def circuit_generator(qarg_gen): # Prepare the initial state for particular problem instance, the default is the \ket{0} state diff --git a/src/qrisp/operators/qubit/__init__.py b/src/qrisp/operators/qubit/__init__.py index 3c956d1a..df4f7e94 100644 --- a/src/qrisp/operators/qubit/__init__.py +++ b/src/qrisp/operators/qubit/__init__.py @@ -20,5 +20,5 @@ from qrisp.operators.qubit.operator_factors import * from qrisp.operators.qubit.bound_qubit_operator import * -from qrisp.operators.qubit.pauli_measurement import * +#from qrisp.operators.qubit.pauli_measurement import * from qrisp.operators.qubit.operator_factors import * diff --git a/src/qrisp/operators/qubit/bound_qubit_operator.py b/src/qrisp/operators/qubit/bound_qubit_operator.py index fe3d37e8..093dffb1 100644 --- a/src/qrisp/operators/qubit/bound_qubit_operator.py +++ b/src/qrisp/operators/qubit/bound_qubit_operator.py @@ -17,7 +17,6 @@ """ from qrisp.operators.hamiltonian import Hamiltonian from qrisp.operators.qubit.bound_qubit_term import BoundQubitTerm -from qrisp.operators.qubit.pauli_measurement import PauliMeasurement from qrisp.operators.qubit.measurement import get_measurement from qrisp import h, sx, IterationEnvironment, conjugate, merge diff --git a/src/qrisp/operators/qubit/measurement.py b/src/qrisp/operators/qubit/measurement.py index bf315f20..af01f1ed 100644 --- a/src/qrisp/operators/qubit/measurement.py +++ b/src/qrisp/operators/qubit/measurement.py @@ -131,40 +131,53 @@ def get_measurement( qc = combine_single_qubit_gates(qc) qc = qc.transpile() - + + hamiltonian = hamiltonian.hermitize() + hamiltonian = hamiltonian.eliminate_ladder_conjugates() + from qrisp.misc import get_measurement_from_qc - - if _measurement is None: - pauli_measurement = hamiltonian.pauli_measurement() - else: - pauli_measurement = _measurement - - meas_circs = pauli_measurement.circuits - meas_ops = pauli_measurement.operators_int - meas_coeffs = pauli_measurement.coefficients - meas_shots = pauli_measurement.shots - - meas_shots = [round(x/precision**2) for x in meas_shots] - tot_shots = sum(x for x in meas_shots) + results = [] + + groups = hamiltonian.commuting_qw_groups() + + # Collect standard deviation + stds = [] + for group in groups: + stds.append(np.sqrt(group.to_pauli().hermitize().get_operator_variance())) + + N = sum(stds) + shots_list = [int(N*s/precision**2) for s in stds] + tot_shots = sum(x for x in shots_list) + if tot_shots>shots: - #meas_shots = [round(x*shots/tot_shots) for x in meas_shots] raise Warning("Warning: The total amount of shots required: " + str(tot_shots) +" for the target precision: " + str(precision) + " exceeds the allowed maxium amount of shots. Decrease the precision or increase the maxium amount of shots.") - results = [] - - for index,circuit in enumerate(meas_circs): - + + meas_ops = [] + meas_coeffs = [] + + for group in groups: + + conjugation_circuit = group.get_conjugation_circuit() + curr = qc.copy() - qubits = [qarg[i] for i in range(circuit.num_qubits())] - curr.append(circuit.to_gate(), qubits) + qubits = [qarg[i] for i in range(conjugation_circuit.num_qubits())] + curr.append(conjugation_circuit.to_gate(), qubits) - # print(curr.transpile()) - res = get_measurement_from_qc(curr.transpile(), list(qarg), backend, meas_shots[index]) + res = get_measurement_from_qc(curr.transpile(), list(qarg), backend, shots_list.pop(0)) results.append(res) - - expectation = evaluate_expectation(results, meas_ops, meas_coeffs) - #expectation = evaluate_expectation_numba(results, meas_ops, meas_coeffs, meas_qubits) - return expectation + + + temp_meas_ops = [] + temp_coeff = [] + for term, coeff in group.terms_dict.items(): + temp_meas_ops.append(term.serialize()) + temp_coeff.append(coeff) + + meas_coeffs.append(temp_coeff) + meas_ops.append(temp_meas_ops) + + return evaluate_expectation(results, meas_ops, meas_coeffs) # # Evaluate expectation @@ -236,7 +249,7 @@ def evaluate_expectation(results, operators, coefficients): for index1,ops in enumerate(operators): for index2,op in enumerate(ops): for outcome,probability in results[index1].items(): - expectation += probability*evaluate_observable(op,outcome)*coefficients[index1][index2] + expectation += probability*evaluate_observable(op,outcome)*np.real(coefficients[index1][index2]) return expectation diff --git a/src/qrisp/operators/qubit/pauli_measurement.py b/src/qrisp/operators/qubit/pauli_measurement_.py similarity index 100% rename from src/qrisp/operators/qubit/pauli_measurement.py rename to src/qrisp/operators/qubit/pauli_measurement_.py diff --git a/src/qrisp/operators/qubit/qubit_operator.py b/src/qrisp/operators/qubit/qubit_operator.py index adc92dd7..3b0ea30d 100644 --- a/src/qrisp/operators/qubit/qubit_operator.py +++ b/src/qrisp/operators/qubit/qubit_operator.py @@ -18,11 +18,11 @@ from qrisp.operators.hamiltonian_tools import find_qw_commuting_groups from qrisp.operators.hamiltonian import Hamiltonian from qrisp.operators.qubit.qubit_term import QubitTerm -from qrisp.operators.qubit.pauli_measurement import PauliMeasurement from qrisp.operators.qubit.measurement import get_measurement from qrisp import h, s, x, IterationEnvironment, conjugate, merge import sympy as sp +import numpy as np threshold = 1e-9 @@ -604,12 +604,10 @@ def adjoint(self): >>> print(H.adjoint()) C_0*A_1*Z_2 """ - res = 0 + new_terms_dict = {} for term, coeff in self.terms_dict.items(): - res += coeff*term.adjoint() - if isinstance(res, (float, int)): - return QubitOperator({QubitTerm({}) : res}) - return res + new_terms_dict[term.adjoint()] = np.conjugate(coeff) + return QubitOperator(new_terms_dict) def hermitize(self): """ @@ -625,6 +623,23 @@ def hermitize(self): """ return 0.5*(self + self.adjoint()) + def eliminate_ladder_conjugates(self): + new_terms_dict = {} + for term, coeff in self.terms_dict.items(): + for factor in term.factor_dict.values(): + if factor in ["A", "C"]: + break + else: + new_terms_dict[term] = coeff + continue + + if term.adjoint() in new_terms_dict: + new_terms_dict[term.adjoint()] += coeff + else: + new_terms_dict[term] = coeff + + return QubitOperator(new_terms_dict) + def ground_state_energy(self): """ @@ -758,10 +773,65 @@ def commuting_qw_groups(self, show_bases=False): # # Measurement settings and measurement # + + def get_conjugation_circuit(self): + + from qrisp import QuantumCircuit + n = self.find_minimal_qubit_amount() + qc = QuantumCircuit(n) + + basis_dict = {} + for term, coeff in self.terms_dict.items(): + + factor_dict = term.factor_dict + for j in range(n): + if j not in factor_dict: + continue + + if j in basis_dict: + assert basis_dict[j] == factor_dict[j] + continue + + if factor_dict[j] not in ["X", "Y", "Z"]: + continue + + basis_dict[j] = factor_dict[j] + + if factor_dict[j]=="X": + qc.h(j) + if factor_dict[j]=="Y": + qc.sx(j) + + ladder_operators = [base for base in term.factor_dict.items() if base[1] in ["A", "C"]] + + if len(ladder_operators): + anchor_factor = ladder_operators[-1] + + if anchor_factor[1] == "A": + qc.x(anchor_factor[0]) + + for j in range(len(ladder_operators)-1): + qc.cx(anchor_factor[0], ladder_operators[j][0]) - def pauli_measurement(self): - return PauliMeasurement(self) + if len(ladder_operators): + if anchor_factor[1] == "A": + qc.x(anchor_factor[0]) + + qc.h(anchor_factor[0]) + + return qc + def get_operator_variance(self): + """ + Calculates the optimal distribution and number of shots following https://quantum-journal.org/papers/q-2021-01-20-385/pdf/. + + """ + var = 0 + for term, coeff in self.terms_dict.items(): + if len(term.factor_dict) != 0: + var += abs(coeff)**2 + return var + def get_measurement( self, qarg, diff --git a/src/qrisp/operators/qubit/qubit_term.py b/src/qrisp/operators/qubit/qubit_term.py index 48d8c954..c18a7564 100644 --- a/src/qrisp/operators/qubit/qubit_term.py +++ b/src/qrisp/operators/qubit/qubit_term.py @@ -34,7 +34,8 @@ class QubitTerm: def __init__(self, factor_dict={}): self.factor_dict = factor_dict - self.hash_value = hash(tuple(sorted(factor_dict.items()))) + + self.hash_value = hash(tuple(sorted(factor_dict.items(), key = lambda x : x[0]))) def update(self, update_dict): self.factor_dict.update(update_dict) @@ -112,26 +113,16 @@ def to_pauli(self): return res def adjoint(self): - from qrisp.operators import X, Y, Z, A, C, P0, P1 - res = 1 - + new_factor_dict = {} for i, factor in self.factor_dict.items(): - if factor == "X": - res *= X(i) - elif factor == "Y": - res *= Y(i) - elif factor == "Z": - res *= Z(i) + if factor in ["X", "Y", "Z", "P0", "P1"]: + new_factor_dict[i] = factor elif factor == "A": - res *= C(i) + new_factor_dict[i] = "C" elif factor == "C": - res *= A(i) - elif factor == "P0": - res *= P0(i) - elif factor == "P1": - res *= P1(i) + new_factor_dict[i] = "A" - return res + return QubitTerm(new_factor_dict) # # Simulation diff --git a/tests/test_VQE.py b/tests/operator_tests/test_VQE.py similarity index 100% rename from tests/test_VQE.py rename to tests/operator_tests/test_VQE.py diff --git a/tests/test_VQE_electronic_structure.py b/tests/operator_tests/test_VQE_electronic_structure.py similarity index 100% rename from tests/test_VQE_electronic_structure.py rename to tests/operator_tests/test_VQE_electronic_structure.py diff --git a/tests/test_VQE_heisenberg.py b/tests/operator_tests/test_VQE_heisenberg.py similarity index 100% rename from tests/test_VQE_heisenberg.py rename to tests/operator_tests/test_VQE_heisenberg.py diff --git a/tests/hamiltonian_tests/test_fermionic_hamiltonian_simulation.py b/tests/operator_tests/test_fermionic_hamiltonian_simulation.py similarity index 100% rename from tests/hamiltonian_tests/test_fermionic_hamiltonian_simulation.py rename to tests/operator_tests/test_fermionic_hamiltonian_simulation.py diff --git a/tests/hamiltonian_tests/test_fermionic_term.py b/tests/operator_tests/test_fermionic_term.py similarity index 100% rename from tests/hamiltonian_tests/test_fermionic_term.py rename to tests/operator_tests/test_fermionic_term.py diff --git a/tests/hamiltonian_tests/test_fermionic_to_qubit.py b/tests/operator_tests/test_fermionic_to_qubit.py similarity index 100% rename from tests/hamiltonian_tests/test_fermionic_to_qubit.py rename to tests/operator_tests/test_fermionic_to_qubit.py diff --git a/tests/test_hamiltonian.py b/tests/operator_tests/test_hamiltonian.py similarity index 100% rename from tests/test_hamiltonian.py rename to tests/operator_tests/test_hamiltonian.py diff --git a/tests/hamiltonian_tests/test_measurement_method.py b/tests/operator_tests/test_measurement_method.py similarity index 93% rename from tests/hamiltonian_tests/test_measurement_method.py rename to tests/operator_tests/test_measurement_method.py index a00c8949..975ae418 100644 --- a/tests/hamiltonian_tests/test_measurement_method.py +++ b/tests/operator_tests/test_measurement_method.py @@ -33,7 +33,7 @@ def testing_helper(qv): continue print(H) - assert abs(H.get_measurement(qv, precision = 0.001) - H.to_pauli().get_measurement(qv, precision = 0.001)) < 1E-2 + assert abs(H.get_measurement(qv, precision = 0.0001, shots = int(1E8)) - H.to_pauli().get_measurement(qv, precision = 0.0001, shots = int(1E8))) < 1E-2 qv = QuantumVariable(4) diff --git a/tests/hamiltonian_tests/test_qubit_hamiltonian_commutator.py b/tests/operator_tests/test_qubit_hamiltonian_commutator.py similarity index 100% rename from tests/hamiltonian_tests/test_qubit_hamiltonian_commutator.py rename to tests/operator_tests/test_qubit_hamiltonian_commutator.py diff --git a/tests/hamiltonian_tests/test_qubit_hamiltonian_simulation.py b/tests/operator_tests/test_qubit_hamiltonian_simulation.py similarity index 100% rename from tests/hamiltonian_tests/test_qubit_hamiltonian_simulation.py rename to tests/operator_tests/test_qubit_hamiltonian_simulation.py diff --git a/tests/hamiltonian_tests/test_to_pauli_converter.py b/tests/operator_tests/test_to_pauli_converter.py similarity index 100% rename from tests/hamiltonian_tests/test_to_pauli_converter.py rename to tests/operator_tests/test_to_pauli_converter.py From 5e9a305c2471286142a2b7da45cf2492f76a9a7d Mon Sep 17 00:00:00 2001 From: positr0nium Date: Tue, 12 Nov 2024 17:26:42 +0100 Subject: [PATCH 07/12] added some comments to QubitOperator.get_conjugation_circuit --- src/qrisp/operators/qubit/measurement.py | 1 - src/qrisp/operators/qubit/qubit_operator.py | 97 +++++++++++++++++++-- 2 files changed, 92 insertions(+), 6 deletions(-) diff --git a/src/qrisp/operators/qubit/measurement.py b/src/qrisp/operators/qubit/measurement.py index af01f1ed..aff4fe0f 100644 --- a/src/qrisp/operators/qubit/measurement.py +++ b/src/qrisp/operators/qubit/measurement.py @@ -167,7 +167,6 @@ def get_measurement( res = get_measurement_from_qc(curr.transpile(), list(qarg), backend, shots_list.pop(0)) results.append(res) - temp_meas_ops = [] temp_coeff = [] for term, coeff in group.terms_dict.items(): diff --git a/src/qrisp/operators/qubit/qubit_operator.py b/src/qrisp/operators/qubit/qubit_operator.py index 3b0ea30d..500c602c 100644 --- a/src/qrisp/operators/qubit/qubit_operator.py +++ b/src/qrisp/operators/qubit/qubit_operator.py @@ -775,48 +775,135 @@ def commuting_qw_groups(self, show_bases=False): # def get_conjugation_circuit(self): + # This method returns a QuantumCircuit that should be applied + # before a measurement of self is peformed. + # The method assumes that all terms within this Operator commute qubit- + # wise. For instance, if an X operator is supposed to be measured, + # the conjugation circuit will contain an H gate at that point, + # because the X operator can be measured by measuring the Z Operator + # in the H-transformed basis. + # For the ladder operators, the conjugation circuit not this straight- + # forward. To understand how we measure the ladder operators, consider + # the operator + + # H = (A(0)*A(1)*A(2) + h.c.) + # = (|000><111| + |111><000|) + + # The considerations from Selingers Paper https://arxiv.org/abs/2310.12256 + # motivate that H can be expressed as a conjugation of the following form. + + # H = U^dg (|110><110| - |111><111|) U + + # This is because + + # exp(i*t*H) = U^dg MCRZ(i*t) U + # = U^dg exp(i*t*(|110><110| - |111><111|)) U + + # We use this insight because the Operator + # |111><111| - |110><110| = |11><11| (x) (|0><0| - |1><1|) + # = |11><11| (x) Z + # can be measured via postprocessing. + + # The postprocessing to do is essentially measuring the last qubit as + # a regular Z operator and only add the result to the expectation value + # if the first two qubits are measured to be in the |1> state. + # If they are in any other state nothing should be added. + + # From this we can also conclude how the conjugation circuit needs to + # look like: Essentially like the conjugation circuit from the paper. + + # For our example above (when simulated) gives: + + # ┌───┐ ┌───┐ + # qv_0.0: ─────┤ X ├────────────■─────────────────────────■─────────────────┤ X ├───── + # └─┬─┘┌───┐ │ │ ┌───┐└─┬─┘ + # qv_0.1: ───────┼──┤ X ├───────■─────────────────────────■────────────┤ X ├──┼─────── + # ┌───┐ │ └─┬─┘┌───┐ │ ┌───┐┌──────────────┐ │ ┌───┐┌───┐└─┬─┘ │ ┌───┐ + # qv_0.2: ┤ X ├──■────■──┤ X ├──┼──┤ H ├┤ Rz(-1.0*phi) ├──┼──┤ H ├┤ X ├──■────■──┤ X ├ + # └───┘ └───┘┌─┴─┐└───┘└──────┬───────┘┌─┴─┐└───┘└───┘ └───┘ + # hs_anc.0: ────────────────────┤ X ├────────────■────────┤ X ├───────────────────────── + # └───┘ └───┘ + + # Where the construction of the MCRZ gate is is encoded into the Toffolis + # and the controlled RZ-Gate. + + # The conjugation circuit therefore needs to look like this: + + # ┌───┐ + # qb_90: ─────┤ X ├─────────────── + # └─┬─┘┌───┐ + # qb_91: ───────┼──┤ X ├────────── + # ┌───┐ │ └─┬─┘┌───┐┌───┐ + # qb_92: ┤ X ├──■────■──┤ X ├┤ H ├ + # └───┘ └───┘└───┘ + + # To learn more about how the post-processing is implemented check the + # comments of QubitTerm.serialize + + # =============== + + # Create a QuantumCircuit that contains the conjugation from qrisp import QuantumCircuit n = self.find_minimal_qubit_amount() qc = QuantumCircuit(n) + # We track which qubit is in which basis to raise an error if a + # violation with the requirement of qubit wise commutativity is detected. basis_dict = {} + + # We iterate through the terms and apply the appropriate basis transformation for term, coeff in self.terms_dict.items(): factor_dict = term.factor_dict for j in range(n): + + # If there is no entry in the factor dict, this corresponds to + # identity => no basis change required. if j not in factor_dict: continue + # If j is already in the basis dict, we assert that the bases agree + # (otherwise there is a violation of qubit-wise commutativity) if j in basis_dict: assert basis_dict[j] == factor_dict[j] continue + # We treat ladder operators in the next section if factor_dict[j] not in ["X", "Y", "Z"]: continue + # Update the basis dict basis_dict[j] = factor_dict[j] + # Append the appropriate basis-change gate if factor_dict[j]=="X": qc.h(j) if factor_dict[j]=="Y": qc.sx(j) - + + # Next we treat the ladder operators ladder_operators = [base for base in term.factor_dict.items() if base[1] in ["A", "C"]] if len(ladder_operators): + + # The anchor factor is the "last" ladder operator. + # This is the qubit where the H gate will be executed. anchor_factor = ladder_operators[-1] + # Flip the anchor qubit if the ladder operator is an annihilator if anchor_factor[1] == "A": qc.x(anchor_factor[0]) - - for j in range(len(ladder_operators)-1): - qc.cx(anchor_factor[0], ladder_operators[j][0]) + + # Perform the cnot gates + for j in range(len(ladder_operators)-1): + qc.cx(anchor_factor[0], ladder_operators[j][0]) - if len(ladder_operators): + # Flip the anchor qubit back if anchor_factor[1] == "A": qc.x(anchor_factor[0]) + # Execute the H-gate qc.h(anchor_factor[0]) return qc From 188c43d52cc2542c5917e0c93ce6f37aa7e4f2e0 Mon Sep 17 00:00:00 2001 From: positr0nium Date: Tue, 12 Nov 2024 17:40:59 +0100 Subject: [PATCH 08/12] wrote some comments for QubitTerm.serialize --- src/qrisp/operators/qubit/measurement.py | 8 +------ src/qrisp/operators/qubit/qubit_term.py | 29 ++++++++++++++++++++++-- 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/src/qrisp/operators/qubit/measurement.py b/src/qrisp/operators/qubit/measurement.py index aff4fe0f..81bce7bc 100644 --- a/src/qrisp/operators/qubit/measurement.py +++ b/src/qrisp/operators/qubit/measurement.py @@ -212,14 +212,8 @@ def evaluate_observable(observable: tuple, x: int): sign_flip = bin(z_int & x).count('1') - from qrisp import bin_rep temp = (x ^ AND_ctrl_state) - # if AND_bits != 0: - # print(bin_rep(AND_ctrl_state, 3)) - # print(bin_rep(AND_bits, 3)) - # print(bin_rep(temp, 3)) - # print(bin_rep(z_int, 3)) - # print("=====") + if contains_ladder: prefactor = 0.5 else: diff --git a/src/qrisp/operators/qubit/qubit_term.py b/src/qrisp/operators/qubit/qubit_term.py index c18a7564..7cd317e8 100644 --- a/src/qrisp/operators/qubit/qubit_term.py +++ b/src/qrisp/operators/qubit/qubit_term.py @@ -54,6 +54,25 @@ def is_identity(self): return len(self.factor_dict)==0 def serialize(self): + # This function serializes the QubitTerm in a way that facilitates the + # measurement post-processing. To learn more details about the reasoning + # behind this Serialisation check QubitOperator.conjugation_circuit + + # The idea here is to serialize the operator via 3 integers. + # These integers specify how the energy of a measurement sample should + # be computed. + # They are processed in the function "evaluate_observable". + + # 1. The Z-int: The binary representation of this integer has a 1 at + # every digit, where there is a Pauli term in self. + + # 2. The and_int: The binary representation has a 1 at every digit, which + # should participate in an AND evaluation. If the evaluation of the AND + # does not return True, the energy of this measurement is 0. + + # 3. The ctrl_int: Thi binary representation has a 1 at every digit, + # which should be flipped before evaluating the AND value. + z_int = 0 and_int = 0 ctrl_int = 0 @@ -62,6 +81,7 @@ def serialize(self): for i in factor_dict.keys(): bit = (1< Date: Tue, 12 Nov 2024 17:48:10 +0100 Subject: [PATCH 09/12] added some comments to evaluate observable --- src/qrisp/operators/qubit/measurement.py | 45 +++++++++--------------- 1 file changed, 16 insertions(+), 29 deletions(-) diff --git a/src/qrisp/operators/qubit/measurement.py b/src/qrisp/operators/qubit/measurement.py index 81bce7bc..ce837ada 100644 --- a/src/qrisp/operators/qubit/measurement.py +++ b/src/qrisp/operators/qubit/measurement.py @@ -184,52 +184,39 @@ def get_measurement( def evaluate_observable(observable: tuple, x: int): - """ - This method evaluates an observable that is a tensor product of Pauli-:math:`Z` operators - with respect to a measurement outcome. - - A Pauli operator of the form :math:`\prod_{i\in I}Z_i`, for some finite set of indices :math:`I\subset \mathbb N`, - is identified with an integer: - We identify the Pauli operator with the binary string that has ones at positions :math:`i\in I` - and zeros otherwise, and then convert this binary string to an integer. - - Parameters - ---------- - - observable : int - The observable represented as integer. - x : int - The measurement outcome represented as integer. - - Returns - ------- - int - The value of the observable with respect to the measurement outcome. - - """ + # This function evaluates how to compute the energy of a measurement sample x. + # Since we are also considering ladder operators, this energy can either be + # 0, -1 or 1. For more details check out the comments of QubitOperator.get_conjugation_circuit - z_int, AND_bits, AND_ctrl_state, contains_ladder = observable + # The observable is given as tuple, containing for integers and a boolean. + # To understand the meaning of these integers check QubitTerm.serialize. + # Unwrap the tuple + z_int, AND_bits, AND_ctrl_state, contains_ladder = observable + + # Compute whether the sign should be sign flipped based on the Z operators sign_flip = bin(z_int & x).count('1') - temp = (x ^ AND_ctrl_state) - + # If the obervable contains ladder operators we need to half the energy? if contains_ladder: prefactor = 0.5 else: prefactor = 1 + # If there are no and bits, we return the result if AND_bits == 0: return prefactor*(-1)**sign_flip + + # Otherwise we apply the AND_ctrl_state to flip the appropriate bits. + corrected_x = (x ^ AND_ctrl_state) - if temp & AND_bits == 0: + # If all bits are in the 0 state the AND is true. + if corrected_x & AND_bits == 0: return prefactor*(-1)**sign_flip else: return 0 - return (-1)**(sign_flip) - def evaluate_expectation(results, operators, coefficients): """ From 09a0e783caa8e1cd5bc47766e234fae8c4a5b846 Mon Sep 17 00:00:00 2001 From: positr0nium Date: Tue, 12 Nov 2024 17:54:28 +0100 Subject: [PATCH 10/12] small improvement of previous comments --- src/qrisp/operators/qubit/measurement.py | 3 ++- src/qrisp/operators/qubit/qubit_operator.py | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/qrisp/operators/qubit/measurement.py b/src/qrisp/operators/qubit/measurement.py index ce837ada..e26a0b8b 100644 --- a/src/qrisp/operators/qubit/measurement.py +++ b/src/qrisp/operators/qubit/measurement.py @@ -197,7 +197,8 @@ def evaluate_observable(observable: tuple, x: int): # Compute whether the sign should be sign flipped based on the Z operators sign_flip = bin(z_int & x).count('1') - # If the obervable contains ladder operators we need to half the energy? + # If there is a ladder operator in the term, we need to half the energy + # because we want to measure (|110><110| - |111><111|)/2 if contains_ladder: prefactor = 0.5 else: diff --git a/src/qrisp/operators/qubit/qubit_operator.py b/src/qrisp/operators/qubit/qubit_operator.py index 500c602c..5648f2ed 100644 --- a/src/qrisp/operators/qubit/qubit_operator.py +++ b/src/qrisp/operators/qubit/qubit_operator.py @@ -793,12 +793,12 @@ def get_conjugation_circuit(self): # The considerations from Selingers Paper https://arxiv.org/abs/2310.12256 # motivate that H can be expressed as a conjugation of the following form. - # H = U^dg (|110><110| - |111><111|) U + # H = U^dg (|110><110| - |111><111|)/2 U # This is because # exp(i*t*H) = U^dg MCRZ(i*t) U - # = U^dg exp(i*t*(|110><110| - |111><111|)) U + # = U^dg exp(i*t*(|110><110| - |111><111|)/2) U # We use this insight because the Operator # |111><111| - |110><110| = |11><11| (x) (|0><0| - |1><1|) From bf3fca91b08cf27fa7db8664482e9545d050dd10 Mon Sep 17 00:00:00 2001 From: positr0nium Date: Tue, 12 Nov 2024 18:46:24 +0100 Subject: [PATCH 11/12] reestablished numba based observable evaluation --- src/qrisp/operators/qubit/measurement.py | 162 ++++++++++------------- src/qrisp/operators/qubit/qubit_term.py | 2 +- 2 files changed, 70 insertions(+), 94 deletions(-) diff --git a/src/qrisp/operators/qubit/measurement.py b/src/qrisp/operators/qubit/measurement.py index e26a0b8b..75056182 100644 --- a/src/qrisp/operators/qubit/measurement.py +++ b/src/qrisp/operators/qubit/measurement.py @@ -176,13 +176,64 @@ def get_measurement( meas_coeffs.append(temp_coeff) meas_ops.append(temp_meas_ops) - return evaluate_expectation(results, meas_ops, meas_coeffs) + + samples = create_padded_array([list(res.keys()) for res in results]).astype(np.int64) + probs = create_padded_array([list(res.values()) for res in results]) + meas_ops = np.array(meas_ops, dtype = np.int64) + meas_coeffs = np.array(meas_coeffs) + + return evaluate_expectation_jitted(samples, probs, meas_ops, meas_coeffs) + + +def create_padded_array(list_of_lists): + """ + Create a padded numpy array from a list of lists with varying lengths. + + Parameters: + list_of_lists (list): A list of lists with potentially different lengths. + + Returns: + numpy.ndarray: A 2D numpy array with padded rows. + """ + # Find the maximum length of any list in the input + max_length = max(len(lst) for lst in list_of_lists) + + # Create a padded list of lists + padded_lists = [ + lst + [0] * (max_length - len(lst)) + for lst in list_of_lists + ] + + # Convert to numpy array + return np.array(padded_lists) + # # Evaluate expectation # + +def evaluate_expectation(samples, probs, operators, coefficients): + """ + Evaluate the expectation. + """ + # print(results) + # print(operators) + # print(coefficients) + # raise + + expectation = 0 + + for index1,ops in enumerate(operators): + for index2,op in enumerate(ops): + for i in range(len(samples[index1])): + outcome,probability = samples[index1, i], probs[index1, i] + expectation += probability*evaluate_observable(op,outcome)*np.real(coefficients[index1][index2]) + + return expectation + + def evaluate_observable(observable: tuple, x: int): # This function evaluates how to compute the energy of a measurement sample x. # Since we are also considering ladder operators, this energy can either be @@ -195,7 +246,11 @@ def evaluate_observable(observable: tuple, x: int): z_int, AND_bits, AND_ctrl_state, contains_ladder = observable # Compute whether the sign should be sign flipped based on the Z operators - sign_flip = bin(z_int & x).count('1') + sign_flip_int = z_int & x + sign_flip = 0 + while sign_flip_int: + sign_flip += sign_flip_int & 1 + sign_flip_int >>= 1 # If there is a ladder operator in the term, we need to half the energy # because we want to measure (|110><110| - |111><111|)/2 @@ -217,24 +272,30 @@ def evaluate_observable(observable: tuple, x: int): else: return 0 - -def evaluate_expectation(results, operators, coefficients): +evaluate_observable_jitted = njit(cache = True)(evaluate_observable) + +@njit(cache = True) +def evaluate_expectation_jitted(samples, probs, operators, coefficients): """ Evaluate the expectation. """ + # print(results) + # print(operators) + # print(coefficients) + # raise expectation = 0 for index1,ops in enumerate(operators): for index2,op in enumerate(ops): - for outcome,probability in results[index1].items(): - expectation += probability*evaluate_observable(op,outcome)*np.real(coefficients[index1][index2]) + for i in range(len(samples[index1])): + outcome,probability = samples[index1, i], probs[index1, i] + expectation += probability*evaluate_observable_jitted(op,outcome)*np.real(coefficients[index1][index2]) return expectation - # # Numba accelearation # @@ -277,89 +338,4 @@ def partition(values, num_qubits): return [np.array(partition[0], dtype=np.uint64)] else: return [np.array(part, dtype=np.uint64) for part in partition[1:]] - - -@njit(cache = True) -def evaluate_observable_jitted(observable, x): - """ - - """ - - value = observable & x - count = 0 - while value: - count += value & 1 - value >>= 1 - return 1 if count % 2 == 0 else -1 - - -@njit(parallel = True, cache = True) -def evaluate_observables_parallel(observables_parts, outcome_parts, probabilities): #M=#observables, N=#measurements - """ - - Parameters - ---------- - observables_parts : list[numpy.array[numpy.unit64]] - The observables. - outcome_parts : list[numpy.array[numpy.unit64]] - The measurement outcomes. - probabilities : numpy.array[numpy.float64] - The measurment probabilities. - - Returns - ------- - - """ - - K = len(observables_parts) # Number of observables PARTS - M = len(observables_parts[0]) # Number of observables - N = len(probabilities) # Number of measurement results - - res = np.zeros(M, dtype=np.float64) - - for j in range(M): - - res_array = np.zeros(N, dtype=np.float64) - for i in prange(N): - temp = 1 - for k in range(K): - temp *= evaluate_observable_jitted(observables_parts[k][j], outcome_parts[k][i]) - - res_array[i] += temp - - res[j] = res_array @ probabilities - - return res - - -def evaluate_expectation_numba(results, operators, coefficients, qubits): - """ - Evaluate the expectation. - - """ - N = len(operators) - - # Partition outcomes in uint64 - #outcomes_parts = [partition(result.keys(), len(meas_qubits[k])) for k, result in enumerate(results)] - outcomes_parts = [partition(list(results[k].keys()), len(qubits[k])) if len(qubits[k])>64 else [np.array(list(results[k].keys()), dtype=np.uint64)] for k in range(N)] - - probabilities = [np.array(list(result.values()), dtype=np.float64) for result in results] - - # Partition observables in uint64 - observables_parts = [partition(observable, len(qubits[k])) if len(qubits[k])>64 else [np.array(observable, dtype=np.uint64)] for k, observable in enumerate(operators)] - - coefficients = [np.array(coeffs, dtype=np.float64) for coeffs in coefficients] - - # Evaluate the observable - expectation = 0 - - - for k in range(N): - - result = evaluate_observables_parallel(observables_parts[k], outcomes_parts[k], probabilities[k]) - #print(result) - #print(coefficients[k]) - - expectation += result @ coefficients[k] - - return expectation \ No newline at end of file + \ No newline at end of file diff --git a/src/qrisp/operators/qubit/qubit_term.py b/src/qrisp/operators/qubit/qubit_term.py index 7cd317e8..194f3519 100644 --- a/src/qrisp/operators/qubit/qubit_term.py +++ b/src/qrisp/operators/qubit/qubit_term.py @@ -112,7 +112,7 @@ def serialize(self): # Returns a tuple, which contains the relevant integers and a boolean # indicating whether this term contains any ladder operators. - return (z_int, and_int, ctrl_int, last_ladder_factor is not None) + return (z_int, and_int, ctrl_int, int(last_ladder_factor is not None)) def to_pauli(self): From 0a3a0ddcf1a07dc2abddaae662df040ade6aca54 Mon Sep 17 00:00:00 2001 From: positr0nium Date: Tue, 12 Nov 2024 19:05:48 +0100 Subject: [PATCH 12/12] fixed a bug in the list padding for numba implementation --- src/qrisp/operators/qubit/measurement.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/src/qrisp/operators/qubit/measurement.py b/src/qrisp/operators/qubit/measurement.py index 75056182..7fcf03f6 100644 --- a/src/qrisp/operators/qubit/measurement.py +++ b/src/qrisp/operators/qubit/measurement.py @@ -179,13 +179,13 @@ def get_measurement( samples = create_padded_array([list(res.keys()) for res in results]).astype(np.int64) probs = create_padded_array([list(res.values()) for res in results]) - meas_ops = np.array(meas_ops, dtype = np.int64) - meas_coeffs = np.array(meas_coeffs) + meas_ops = create_padded_array(meas_ops, use_tuples = True).astype(np.int64) + meas_coeffs = create_padded_array(meas_coeffs) return evaluate_expectation_jitted(samples, probs, meas_ops, meas_coeffs) -def create_padded_array(list_of_lists): +def create_padded_array(list_of_lists, use_tuples = False): """ Create a padded numpy array from a list of lists with varying lengths. @@ -199,10 +199,17 @@ def create_padded_array(list_of_lists): max_length = max(len(lst) for lst in list_of_lists) # Create a padded list of lists - padded_lists = [ - lst + [0] * (max_length - len(lst)) - for lst in list_of_lists - ] + if not use_tuples: + padded_lists = [ + lst + [0] * (max_length - len(lst)) + for lst in list_of_lists + ] + else: + padded_lists = [ + lst + [(0,0,0,0)] * (max_length - len(lst)) + for lst in list_of_lists + ] + # Convert to numpy array return np.array(padded_lists)