From b2fc104519e65152887bc5603396a933cc93c2c7 Mon Sep 17 00:00:00 2001 From: Nick Rubin Date: Tue, 17 Jul 2018 15:55:44 +1000 Subject: [PATCH 01/10] WIP stabilizer --- referenceqvm/api.py | 2 + referenceqvm/gates.py | 10 +++ referenceqvm/qvm_stabilizer.py | 112 +++++++++++++++++++++++++++++++++ 3 files changed, 124 insertions(+) create mode 100644 referenceqvm/qvm_stabilizer.py diff --git a/referenceqvm/api.py b/referenceqvm/api.py index 422770c..5d6d4db 100644 --- a/referenceqvm/api.py +++ b/referenceqvm/api.py @@ -35,6 +35,8 @@ def QVMConnection(type_trans='wavefunction', elif type_trans == 'density': qvm = QVM_Density(gate_set=gate_set, noise_model=noise_model) + elif type_trans == 'stabilizer': + qvm = QVM_Stabilizer(gate_set=stabilizer_gate_set) else: raise TypeError("{} is not a valid QVM type.".format(type_trans)) diff --git a/referenceqvm/gates.py b/referenceqvm/gates.py index 6cf84be..ca3c833 100644 --- a/referenceqvm/gates.py +++ b/referenceqvm/gates.py @@ -226,6 +226,16 @@ def BARENCO(alpha, phi, theta): 'CZ': CZ } +stabilizer_gate_matrix = { + 'I': I, + 'X': X, + 'Y': Y, + 'Z': Z, + 'H': H, + 'CNOT': CNOT, + 'S': S, + 'CZ': CZ +} # noisy gates def relaxation(p): diff --git a/referenceqvm/qvm_stabilizer.py b/referenceqvm/qvm_stabilizer.py new file mode 100644 index 0000000..cd2d9e4 --- /dev/null +++ b/referenceqvm/qvm_stabilizer.py @@ -0,0 +1,112 @@ +#!/usr/bin/python +############################################################################## +# Copyright 2016-2017 Rigetti Computing +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +############################################################################## +""" +Pure QVM that only executes pyQuil programs containing stabilizer group gates, +and returns the unitary resulting from the program evolution. +""" +from pyquil.quil import Program +from pyquil.quilbase import * + +from referenceqvm.unitary_generator import tensor_gates +from referenceqvm.qam import QAM + + +class QVM_Unitary(QAM): + """ + A P Y T H O N + + Q U A N T U M + + V I R T U A L + + M A C H I N E + + Only pyQuil programs containing pure Gates or DefGate objects are accepted. + The QVM_Unitary kernel applies all the gates, and returns the unitary + corresponding to the input program. + + Note: no classical control flow, measurements allowed. + """ + def __init__(self, num_qubits=None, program=None, program_counter=None, + classical_memory=None, gate_set=None, defgate_set=None, + unitary=None): + """ + Subclassed from QAM this is a pure QVM. + """ + super(QVM_Unitary, self).__init__(num_qubits=num_qubits, program=program, + program_counter=program_counter, + classical_memory=classical_memory, + gate_set=gate_set, + defgate_set=defgate_set) + self.umat = unitary + self.all_inst = False + + def transition(self, instruction): + """ + Implements a transition on the unitary-qvm. + + :param Gate instruction: QuilAction gate to be implemented + """ + if instruction.name in self.gate_set or instruction.name in self.defgate_set: + # get the unitary and evolve the state + unitary = tensor_gates(self.gate_set, self.defgate_set, + instruction, self.num_qubits) + self.umat = unitary.dot(self.umat) + self.program_counter += 1 + else: + raise TypeError("Gate {} is not in the " + "gate set".format(instruction.name)) + + def unitary(self, pyquil_program): + """ + Return the unitary of a pyquil program + + This method initializes a qvm with a gate_set, protoquil program (expressed + as a pyquil program), and then executes the QVM statemachine. + + :param pyquil_program: (pyquil.Program) object containing only protoQuil + instructions. + + :return: a unitary corresponding to the output of the program. + :rtype: np.array + """ + + # load program + self.load_program(pyquil_program) + + # setup unitary + self.umat = np.eye(2 ** self.num_qubits) + + # evolve unitary + self.kernel() + + return self.umat + + def expectation(self, pyquil_program, operator_programs=[Program()]): + """ + Calculate the expectation value given a state prepared. + + :param pyquil_program: (pyquil.Program) object containing only protoQuil + instructions. + :param operator_programs: (optional, list) of PauliTerms. Default is + Identity operator. + + :return: expectation value of the operators. + :rtype: float + """ + # TODO + raise NotImplementedError() From be741eb633b42d325cf5ee8b59ccbd1736e378a4 Mon Sep 17 00:00:00 2001 From: Nick Rubin Date: Sat, 28 Jul 2018 03:06:21 +1000 Subject: [PATCH 02/10] WIP CHP stabilizer sim --- referenceqvm/qvm_stabilizer.py | 480 +++++++++++++++++++--- referenceqvm/state_actions.py | 39 ++ referenceqvm/tests/test_stabilzier_qvm.py | 291 +++++++++++++ referenceqvm/unitary_generator.py | 6 +- 4 files changed, 761 insertions(+), 55 deletions(-) create mode 100644 referenceqvm/state_actions.py create mode 100644 referenceqvm/tests/test_stabilzier_qvm.py diff --git a/referenceqvm/qvm_stabilizer.py b/referenceqvm/qvm_stabilizer.py index cd2d9e4..5e6c39e 100644 --- a/referenceqvm/qvm_stabilizer.py +++ b/referenceqvm/qvm_stabilizer.py @@ -15,17 +15,21 @@ # limitations under the License. ############################################################################## """ -Pure QVM that only executes pyQuil programs containing stabilizer group gates, -and returns the unitary resulting from the program evolution. +Pure QVM that only executes pyQuil programs containing Clifford group generators, +and return the wavefunction or stabilizer """ +import sys +from functools import reduce from pyquil.quil import Program from pyquil.quilbase import * +from pyquil.paulis import PauliTerm, sI, sZ, sX, sY -from referenceqvm.unitary_generator import tensor_gates from referenceqvm.qam import QAM +from referenceqvm.gates import stabilizer_gate_matrix +from referenceqvm.unitary_generator import value_get -class QVM_Unitary(QAM): +class QVM_Stabilizer(QAM): """ A P Y T H O N @@ -35,78 +39,448 @@ class QVM_Unitary(QAM): M A C H I N E - Only pyQuil programs containing pure Gates or DefGate objects are accepted. - The QVM_Unitary kernel applies all the gates, and returns the unitary - corresponding to the input program. + S I M U L A T I N G - Note: no classical control flow, measurements allowed. + S T A B I L I Z E R + + S T A T E S """ def __init__(self, num_qubits=None, program=None, program_counter=None, - classical_memory=None, gate_set=None, defgate_set=None, - unitary=None): + classical_memory=None, gate_set=stabilizer_gate_matrix, + defgate_set=None): """ Subclassed from QAM this is a pure QVM. """ - super(QVM_Unitary, self).__init__(num_qubits=num_qubits, program=program, - program_counter=program_counter, - classical_memory=classical_memory, - gate_set=gate_set, - defgate_set=defgate_set) - self.umat = unitary - self.all_inst = False + super(QVM_Stabilizer, self).__init__(num_qubits=num_qubits, program=program, + program_counter=program_counter, + classical_memory=classical_memory, + gate_set=gate_set, + defgate_set=defgate_set) - def transition(self, instruction): + if num_qubits is None: + self.tableau = None + else: + self.tableau = self._n_qubit_tableau(num_qubits) + + def _n_qubit_tableau(self, num_qubits): """ - Implements a transition on the unitary-qvm. + Construct an empty tableau for a n-qubit system - :param Gate instruction: QuilAction gate to be implemented + :param num_qubits: + :return: """ - if instruction.name in self.gate_set or instruction.name in self.defgate_set: - # get the unitary and evolve the state - unitary = tensor_gates(self.gate_set, self.defgate_set, - instruction, self.num_qubits) - self.umat = unitary.dot(self.umat) - self.program_counter += 1 + tableau = np.zeros((2 * num_qubits, (2 * num_qubits) + 1), + dtype=int) + + # set up initial state |0>^{otimes n} + for ii in range(2 * self.num_qubits): + tableau[ii, ii] = 1 + return tableau + + def run(self, pyquil_program): + """ + Run comman + + :param porgram: + :return: + """ + self.load_program(pyquil_program) + + # set up stabilizers + self.tableau = self._n_qubit_tableau(self.num_qubits) + + def stabilizer_tableau(self): + """ + return the stabilizer part of the tableau + + :return: stabilizer matrix + """ + return self.tableau[self.num_qubits:, :] + + def destabilizer_tableau(self): + """ + Return the destabilizer part of the tableau + + :return: destabilizer matrix + """ + return self.tableau[:self.num_qubits, :] + + def _rowsum(self, h, i): + """ + Implementation of Aaronson-Gottesman rowsum algorithm + + :param Int h: row index 'h' + :param Int i: row index 'i' + :return: + """ + # NOTE: this is left multiplication of P(i). P(i) * P(h) + phase_accumulator = self._rowsum_phase_accumulator(h, i) + # now set the r_{h} + if phase_accumulator % 4 == 0: + self.tableau[h, -1] = 0 + elif phase_accumulator % 4 == 2: + self.tableau[h, -1] = 1 else: - raise TypeError("Gate {} is not in the " - "gate set".format(instruction.name)) + raise ValueError("An impossible value for the phase_accumulator has occurred") + + # now update the rows of the tableau + for j in range(self.num_qubits): + self.tableau[h, j] = self.tableau[i, j] ^ self.tableau[h, j] + self.tableau[h, j + self.num_qubits] = self.tableau[i, j + self.num_qubits] ^ \ + self.tableau[h, j + self.num_qubits] - def unitary(self, pyquil_program): + def _rowsum_phase_accumulator(self, h, i): """ - Return the unitary of a pyquil program + phase accumulator sub algorithm for the rowsum routine - This method initializes a qvm with a gate_set, protoquil program (expressed - as a pyquil program), and then executes the QVM statemachine. + Note: this accumulates the $i$ phase for row_i * row_h NOT row_h, row_i + :param Int h: row index 'h' + :param Int i: row index 'i' + :return: phase mod 4 + """ + phase_accumulator = 0 + for j in range(self.num_qubits): + # self._g_update(x_{hj}, z_{hj}, x_{ij}, z_{ij}) + phase_accumulator += self._g_update(self.tableau[i, j], + self.tableau[i, self.num_qubits + j], + self.tableau[h, j], + self.tableau[h, self.num_qubits + j]) - :param pyquil_program: (pyquil.Program) object containing only protoQuil - instructions. + phase_accumulator += 2 * self.tableau[h, -1] + phase_accumulator += 2 * self.tableau[i, -1] + return phase_accumulator % 4 - :return: a unitary corresponding to the output of the program. - :rtype: np.array + def _g_update(self, x1, z1, x2, z2): """ + function that takes 4 bits and returns the power of $i$ {0, 1, -1} - # load program - self.load_program(pyquil_program) + when the pauli matrices represented by x1z1 and x2z2 are multiplied together + + :param x1: + :param z1: + :param x2: + :param z2: + :return: + """ + # if the first term is identity + if x1 == z1 == 0: + return 0 + + # if the first term is Y + # Y * Z = (1 - 0) = 1j ^ { 1} = 1j + # Y * I = (0 - 0) = 1j ^ { 0} = 1 + # Y * X = (0 - 1) = 1j ^ {-1} = -1j + # Y * Y = (1 - 1) = 1j ^ { 0} = 1 + if x1 == z1 == 1: + return z2 - x2 + + # if the first term is X return z2 * (2 * x2 - 1) + # X * I = (0 * (2 * 0 - 1) = 0 -> 1j^{ 0} = 1 + # X * X = (0 * (2 * 1 - 1) = 0 -> 1j^{ 0} = 1 + # X * Y = (1 * (2 * 1 - 1) = 1 -> 1j^{ 1} = 1j + # X * Z = (1 * (2 * 0 - 1) = -1 -> 1j^{-1} = -1j + if x1 == 1 and z1 == 0: + return z2 * (2 * x2 - 1) + + # if the first term is Z return x2 * (1 - 2 * z2) + # Z * I = (0 * (1 - 2 * 0)) = 0 -> 1j^{ 0} = 1 + # Z * X = (1 * (1 - 2 * 0)) = 1 -> 1j^{ 1} = 1j + # Z * Y = (1 * (1 - 2 * 1)) = -1 -> 1j^{-1} = -1j + # Z * Z = (0 * (1 - 2 * 1)) = 0 -> 1j^{0} = 1 + if x1 == 0 and z1 == 1: + return x2 * (1 - 2 * z2) + + raise ValueError("we were unable to multiply the pauli operators together!") + + def _apply_cnot(self, instruction): + """ + Apply a CNOT to the tableau + + :param instruction: pyquil abstract instruction. Must have + """ + a, b = list(instruction.get_qubits()) # control (a) and target (b) + for i in range(2 * self.num_qubits): + self.tableau[i, -1] = self._cnot_phase_update(i, a, b) + self.tableau[i, b] = self.tableau[i, b] ^ self.tableau[i, a] + self.tableau[i, a + self.num_qubits] = self.tableau[i, a + self.num_qubits] ^ self.tableau[i, b + self.num_qubits] + + def _cnot_phase_update(self, i, c, t): + """ + update r_{i} - # setup unitary - self.umat = np.eye(2 ** self.num_qubits) + :param i: tableau row index + :param c: control qubit index + :param t: target qubit index + :return: 0/1 phase update for r_{i} + """ + return self.tableau[i, -1] ^ (self.tableau[i, c] * self.tableau[i, t + self.num_qubits]) * ( + self.tableau[i, t] ^ self.tableau[i, c + self.num_qubits] ^ 1) - # evolve unitary - self.kernel() + def _apply_hadamard(self, instruction): + """ + Apply a hadamard gate on qubit defined in instruction - return self.umat + :param instruction: + :return: + """ + qubit_label = list(instruction.get_qubits())[0] + for i in range(2 * self.num_qubits): + self.tableau[i, -1] = self.tableau[i, -1] ^ (self.tableau[i, qubit_label] * self.tableau[i, qubit_label + self.num_qubits]) + self.tableau[i, [qubit_label, qubit_label + self.num_qubits]] = self.tableau[i, [qubit_label + self.num_qubits, qubit_label]] - def expectation(self, pyquil_program, operator_programs=[Program()]): + def _apply_phase(self, instruction): """ - Calculate the expectation value given a state prepared. + Apply the phase gate instruction ot the tableau + + :param instruction: + :return: + """ + qubit_label = list(instruction.get_qubits())[0] + for i in range(2 * self.num_qubits): + self.tableau[i, -1] = self.tableau[i, -1] ^ (self.tableau[i, qubit_label] * self.tableau[i, qubit_label + self.num_qubits]) + self.tableau[i, qubit_label + self.num_qubits] = self.tableau[i, qubit_label + self.num_qubits] ^ self.tableau[i, qubit_label] + + def _apply_measurement(self, instruction): + t_qbit = value_get(instruction.qubit) + t_cbit = value_get(instruction.classical_reg) + + print(self.tableau) + check_xpa = False # check if x_pa = 1 for p in {n + 1 ... 2 n} of tableau + if any(self.tableau[self.num_qubits:, t_qbit] == 1): + # find the first one. + xpa_idx = np.where(self.tableau[self.num_qubits:, t_qbit] == 1)[0][0] # take the first index + xpa_idx += self.num_qubits # adjust so we can index into the tableau + for ii in range(2 * self.num_qubits): # loop over each row and call rowsum(ii, xpa_idx) + if ii != xpa_idx and self.tableau[ii, t_qbit] == 1: + self._rowsum(ii, xpa_idx) + self.tableau[xpa_idx - self.num_qubits, :] = self.tableau[xpa_idx, :] + + self.tableau[xpa_idx, :] = np.zeros((1, 2 * self.num_qubits + 1)) + self.tableau[xpa_idx, t_qbit + self.num_qubits] = 1 + # perform the measurement + self.tableau[xpa_idx, -1] = 1 if np.random.random() > 0.5 else 0 + + # set classical results to return + self.classical_memory[t_cbit] = self.tableau[xpa_idx, -1] + + else: + # augment tableaue with a scratch space + self.tableau = np.vstack((self.tableau, np.zeros((1, 2 * self.num_qubits + 1)))) + for ii in range(self.num_qubits): + self._rowsum(2 * self.num_qubits + 1, ii) + + # set the classical bit to be the last element of the scratch row + self.classical_memory[t_cbit] = self.tableau[-1, -1] - :param pyquil_program: (pyquil.Program) object containing only protoQuil - instructions. - :param operator_programs: (optional, list) of PauliTerms. Default is - Identity operator. + # remove the scratch space + self.tableau = self.tableau[:2 * self.num_qubits, :] - :return: expectation value of the operators. - :rtype: float + def transition(self, instruction): + """ + Implements a transition on the generator matrix representing the stabilizers + + :param Gate instruction: QuilAction gate to be implemented """ - # TODO - raise NotImplementedError() + if isinstance(instruction, Measurement): + # mutates classical memory! I should change this... + self._apply_measurement(instruction) + self.program_counter += 1 + + elif isinstance(instruction, Gate): + # apply Gate or DefGate + if Gate.name == 'H': + self._apply_hadamard(instruction) + elif Gate.name == 'S': + self._apply_phase(instruction) + elif Gate.name == 'CNOT': + self._apply_cnot(instruction) + self.program_counter += 1 + + elif isinstance(instruction, Jump): + # unconditional Jump; go directly to Label + self.program_counter = self.find_label(instruction.target) + + elif isinstance(instruction, JumpTarget): + # Label; pass straight over + self.program_counter += 1 + + +def pauli_stabilizer_to_binary_stabilizer(stabilizer_list): + """ + Convert a list of stabilizers represented as PauliTerms to a binary tableau form + + :param List stabilizer_list: list of stabilizers where each element is a PauliTerm + :return: return an integer matrix representing the stabilizers where each row is a + stabilizer. The size of the matrix is n x (2 * n) where n is the maximum + qubit index. + """ + if not all([isinstance(x, PauliTerm) for x in stabilizer_list]): + raise TypeError("At least one element of stabilizer_list is not a PauliTerm") + + max_qubit = max([max(x.get_qubits()) for x in stabilizer_list]) + 1 + stabilizer_tableau = np.zeros((len(stabilizer_list), 2 * max_qubit + 1), dtype=int) + for row_idx, term in enumerate(stabilizer_list): + for i, pterm in term: # iterator for each tensor-product element of the Pauli operator + if pterm == 'X': + stabilizer_tableau[row_idx, i] = 1 + elif pterm == 'Z': + stabilizer_tableau[row_idx, i + max_qubit] = 1 + elif pterm == 'Y': + stabilizer_tableau[row_idx, i] = 1 + stabilizer_tableau[row_idx, i + max_qubit] = 1 + else: + # term is identity + pass + + if not (np.isclose(term.coefficient, -1) or np.isclose(term.coefficient, 1)): + raise ValueError("stabilizers must have a +/- coefficient") + + if int(np.sign(term.coefficient.real)) == 1: + stabilizer_tableau[row_idx, -1] = 0 + elif int(np.sign(term.coefficient.real)) == -1: + stabilizer_tableau[row_idx, -1] = 1 + else: + raise TypeError('unrecognized on pauli term of stabilizer') + + return stabilizer_tableau + + +def binary_stabilizer_to_pauli_stabilizer(stabilizer_tableau): + """ + Convert a stabilizer tableau to a list of PauliTerms + + :param stabilizer_tableau: + :return: + """ + stabilizer_list = [] + num_qubits = (stabilizer_tableau.shape[1] - 1) // 2 + for nn in range(stabilizer_tableau.shape[0]): # iterate through the rows + stabilizer_element = [] + for ii in range(num_qubits): + if stabilizer_tableau[nn, ii] == 1 and stabilizer_tableau[nn, ii + num_qubits] == 0: + stabilizer_element.append(sX(ii)) + elif stabilizer_tableau[nn, ii] == 0 and stabilizer_tableau[nn, ii + num_qubits] == 1: + stabilizer_element.append(sZ(ii)) + elif stabilizer_tableau[nn, ii] == 1 and stabilizer_tableau[nn, ii + num_qubits] == 1: + stabilizer_element.append(sY(ii)) + + stabilizer_term = reduce(lambda x, y: x * y, stabilizer_element) * ((-1) ** stabilizer_tableau[nn, -1]) + stabilizer_list.append(stabilizer_term) + return stabilizer_list + + +def binary_rref(code_matrix): + """ + Convert the binary code_matrix into rref + + :param code_matrix: + :return: + """ + code_matrix = code_matrix.astype(int) + rdim, cdim = code_matrix.shape + for ridx in range(rdim): + # print("start [{},{}] = {}".format(ridx, ridx, + # code_matrix[ridx, ridx])) + # set first element + if not np.isclose(code_matrix[ridx, ridx], 1): + for ii in range(ridx + 1, rdim): + if np.isclose(code_matrix[ii, ridx], 1.0): + # switch rows + # code_matrix[ridx, :], code_matrix[ii, :] = code_matrix[ii, :], code_matrix[ridx, :] + code_matrix[[ridx, ii]] = code_matrix[[ii, ridx]] + print("r_{} <-> r_{}".format(ii + 1, ridx + 1)) + + # print(code_matrix) + # start elimination of all other rows + # print("pivot [{}, {}] = {}".format(ridx, ridx, + # code_matrix[ridx, ridx])) + + rows_to_eliminate = set(range(rdim)) + rows_to_eliminate.remove(ridx) + # print(rows_to_eliminate) + for elim_ridx in list(rows_to_eliminate): + # eliminate by Hadamard product on the rows (XOR element wise) + # if other row element is 1 then eliminate by subtracting ridx row + + # if np.isclose(code_matrix[elim_ridx, ridx], 1.0): + # code_matrix[elim_ridx, :] = code_matrix[elim_ridx, :] - \ + # code_matrix[ridx, :] + if np.isclose(code_matrix[elim_ridx, ridx], 1.0): + code_matrix[elim_ridx, :] = code_matrix[elim_ridx, :] ^ code_matrix[ridx, :] + # print("elim [{}, :]".format(elim_ridx)) + + # print("semi rref at row {}".format(ridx)) + # print(code_matrix) + return code_matrix + + +if __name__ == "__main__": + # practice my clifford transformations + # = |0000...0N> + # N |psi> = N g |psi> = NgNd N|psi> = g2 |psi2> = |psi2> + + from referenceqvm.gates import X, Y, Z, H, CNOT, S, I + from sympy import Matrix + + # HZH = X + test_X = H.dot(Z).dot(H) + assert np.allclose(test_X, X) + + # HXH = Z + test_Z = H.dot(X).dot(H) + assert np.allclose(test_Z, Z) + + # HYH = iHZHHXH = iXZ = -Y + test_nY = H.dot(Y).dot(H) + assert np.allclose(test_nY, -Y) + + # XYX = -Y + assert np.allclose(X.dot(Y).dot(X), -Y) + + # XZX = -Z + assert np.allclose(X.dot(Z).dot(X), -Z) + + # SXS = Y + assert np.allclose(S.dot(X).dot(np.conj(S).T), Y) + + # CNOT(X otimes X)CNOT = X otimes I + # (C, T) for CNOT as written in reference-qvm + assert np.allclose(CNOT.dot(np.kron(I, X)).dot(CNOT), np.kron(I, X)) + assert np.allclose(CNOT.dot(np.kron(X, I)).dot(CNOT), np.kron(X, X)) + assert np.allclose(CNOT.dot(np.kron(X, X)).dot(CNOT), np.kron(X, I)) + + # CNOT(Z otimes Z)CNOT = Z otimes I + assert np.allclose(CNOT.dot(np.kron(Z, Z)).dot(CNOT), np.kron(I, Z)) + assert np.allclose(CNOT.dot(np.kron(Z, I)).dot(CNOT), np.kron(Z, I)) + assert np.allclose(CNOT.dot(np.kron(I, Z)).dot(CNOT), np.kron(Z, Z)) + + Z1 = np.kron(Z, I) + Z2 = np.kron(I, Z) + ZZ = np.kron(Z, Z) + X1 = np.kron(X, I) + X2 = np.kron(I, X) + XX = np.kron(X, X) + YY = np.kron(Y, Y) + + generator_matrix = np.zeros((3, 4)) + generator_matrix[0, 2] = generator_matrix[0, 3] = 1 # ZZ + generator_matrix[1, 0] = generator_matrix[1, 1] = 1 # XX + generator_matrix[2, 0] = generator_matrix[2, 2] = 1 # Y + generator_matrix[2, 1] = generator_matrix[2, 3] = 1 # Y + parity_vector = np.zeros((3, 1)) + parity_vector[-1] = 1 + + # define the amplitudes from the generator set by diagonalizing + # the tableau matrix? + generator_matrix = np.hstack((generator_matrix, parity_vector)) + print(generator_matrix) + + generator_matrix_2 = binary_rref(generator_matrix) + print('\n') + print(generator_matrix_2) + + mat = Matrix(generator_matrix) + mat, pivots = mat.rref() + mat = np.asarray(mat).astype(int) + print(mat) + diff --git a/referenceqvm/state_actions.py b/referenceqvm/state_actions.py new file mode 100644 index 0000000..2e85709 --- /dev/null +++ b/referenceqvm/state_actions.py @@ -0,0 +1,39 @@ +""" +State actions on classical states + +Commonly, in order to recover a state you need to compute the action +of Pauli operators on classical basis states. + +In this module we provide infrastructure to do this for Pauli Operators from +pyquil. + +Given +""" +import numpy as np +from pyquil.paulis import PauliTerm + + +def compute_action(classical_state, pauli_operator, num_qubits): + """ + Compute action of Pauli opertors on a classical state + + :param classical_state: binary repr of a state or an integer. Should be + left most bit (0th position) is the most significant bit + :param num_qubits: + :return: + """ + if not isinstance(pauli_operator, PauliTerm): + raise TypeError("pauli_operator must be a PauliTerm") + + if not isinstance(classical_state, (list, int)): + raise TypeError("classical state must be a list or an integer") + + if isinstance(classical_state, int): + if classical_state < 0: + raise TypeError("classical_state must be a positive integer") + + classical_state = list(map(int, np.binary_repr(classical_state, + width=num_qubits))) + if len(classical_state) != num_qubits: + raise TypeError("classical state not long enough") + diff --git a/referenceqvm/tests/test_stabilzier_qvm.py b/referenceqvm/tests/test_stabilzier_qvm.py new file mode 100644 index 0000000..0a79728 --- /dev/null +++ b/referenceqvm/tests/test_stabilzier_qvm.py @@ -0,0 +1,291 @@ +""" +Testing the QVM stablizer +""" +import sys +import numpy as np +from itertools import product +from functools import reduce +from pyquil.paulis import sI, sX, sY, sZ, PAULI_COEFF, PAULI_OPS +from pyquil.quil import Program +from pyquil.gates import H, S, CNOT +from referenceqvm.qvm_stabilizer import (QVM_Stabilizer, pauli_stabilizer_to_binary_stabilizer, + binary_stabilizer_to_pauli_stabilizer) +from referenceqvm.unitary_generator import value_get + +pauli_subgroup = [sI, sX, sY, sZ] +five_qubit_code_generators = [sX(0) * sZ(1) * sZ(2) * sX(3) * sI(4), + sI(0) * sX(1) * sZ(2) * sZ(3) * sX(4), + sX(0) * sI(1) * sX(2) * sZ(3) * sZ(4), + sZ(0) * sX(1) * sI(2) * sX(3) * sZ(4)] +bell_stabilizer = [sZ(0) * sZ(1), sX(0) * sX(1)] + + +def test_initialization(): + """ + Test if upon initialization the correct size tableau is set up + """ + num_qubits = 4 + qvmstab = QVM_Stabilizer(num_qubits=num_qubits) + assert qvmstab.num_qubits == num_qubits + assert qvmstab.tableau.shape == (2 * num_qubits, 2 * num_qubits + 1) + + initial_tableau = np.hstack((np.eye(2 * num_qubits), np.zeros((2 * num_qubits, 1)))) + assert np.allclose(initial_tableau, qvmstab.tableau) + + num_qubits = 6 + qvmstab = QVM_Stabilizer(num_qubits=num_qubits) + assert qvmstab.num_qubits == num_qubits + assert qvmstab.tableau.shape == (2 * num_qubits, 2 * num_qubits + 1) + + initial_tableau = np.hstack((np.eye(2 * num_qubits), np.zeros((2 * num_qubits, 1)))) + assert np.allclose(initial_tableau, qvmstab.tableau) + + +def test_stabilizer_to_matrix_conversion(): + # bitflip code + stabilizer_matrix = pauli_stabilizer_to_binary_stabilizer(bell_stabilizer) + true_stabilizer_matrix = np.array([[0, 0, 1, 1, 0], + [1, 1, 0, 0, 0]]) + assert np.allclose(true_stabilizer_matrix, stabilizer_matrix) + + test_stabilizer_list = binary_stabilizer_to_pauli_stabilizer(true_stabilizer_matrix) + for idx, term in enumerate(test_stabilizer_list): + assert term == bell_stabilizer[idx] + + # given some codes convert them to code matrices + stabilizer_matrix = pauli_stabilizer_to_binary_stabilizer(five_qubit_code_generators) + true_stabilizer_matrix = np.array([[1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0], + [0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0], + [1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0]]) + assert np.allclose(true_stabilizer_matrix, stabilizer_matrix) + + test_stabilizer_list = binary_stabilizer_to_pauli_stabilizer(true_stabilizer_matrix) + for idx, term in enumerate(test_stabilizer_list): + assert term == five_qubit_code_generators[idx] + + +def test_row_sum_sub_algorithm_g_test(): + """ + Test the behavior of the subroutine of the rowsum routine + + g(x1, z1, x2, z2) is a function that returns the exponent that $i$ is raised (0, -1, 1) + when pauli matrices represented by x1z1 and x2z2 are multiplied together. + + Logic table: + I1 = {x1 = 0, z1 = 0} + I2 = {x2 = 0, z2 = 0} + g(I1, I2) = 0 + + I1 = {x1 = 0, z1 = 0} + X2 = {x2 = 1, z2 = 0} + g(I1, X2) = 0 + + I1 = {x1 = 0, z1 = 0} + Y2 = {x2 = 1, z2 = 1} + g(I1, X2) = 0 + + I1 = {x1 = 0, z1 = 0} + Z2 = {x2 = 0, z2 = 1} + g(I1, X2) = 0 + + --------------------- + + X1 = {x1 = 1, z1 = 0} + I2 = {x2 = 0, z2 = 0} + g(I1, I2) = 0 + + X1 = {x1 = 1, z1 = 0} + X2 = {x2 = 1, z2 = 0} + g(I1, X2) = 0 + + X1 = {x1 = 1, z1 = 0} + Y2 = {x2 = 1, z2 = 1} + g(I1, X2) = 1 + + X1 = {x1 = 1, z1 = 0} + Z2 = {x2 = 0, z2 = 1} + g(I1, X2) = -1 + + --------------------- + etc... + """ + num_qubits = 4 + qvmstab = QVM_Stabilizer(num_qubits=num_qubits) + + # loop over the one-qubit pauli group and make sure the power of $i$ is correct + # this test would be better if we stored levi-civta + pauli_map = {'X': (1, 0), 'Y': (1, 1), 'Z': (0, 1), 'I': (0, 0)} + for pi, pj in product(PAULI_OPS, repeat=2): + i_coeff = qvmstab._g_update(pauli_map[pi][0], pauli_map[pi][1], pauli_map[pj][0], pauli_map[pj][1]) + # print(pi, pj, "".join([x for x in [pi, pj]]), PAULI_COEFF["".join([x for x in [pi, pj]])], 1j**i_coeff) + assert np.isclose(PAULI_COEFF["".join([x for x in [pi, pj]])], 1j**i_coeff) + + +def test_rowsum_phase_accumulator(): + """ + Test the accuracy of the phase accumulator. This subroutine keeps track of the power that $i$ is raised to + when multiplying two PauliTerms together. PauliTerms are now composed of multi single-qubit Pauli's. + + The formula is phase_accumulator(h, i) = 2 * rh + 2 * ri + \sum_{j}^{n}g(x_{ij}, z_{ij}, x_{hj}, z_{hj} + + The returned value is presented mod 4. + + notice that since r_h indicates +1 or -1 this corresponds to i^{0} or i^{2}. Given r_{h/i} takes on {0, 1} then + we need to multiply by 2 get the correct power of $i$ for the Pauli Term to account for multiplication. + + In order to test we will generate random elements of P_{n} and then multiply them together keeping track of the + $i$ power. We will then compare this to the phase_accumulator subroutine. + + We have to load in the stabilizer into an empty qvmstab object because the subroutine needs a tableau as a reference + """ + num_qubits = 2 + pauli_terms = [sX(0) * sX(1), sZ(0) * sZ(1)] + stab_mat = pauli_stabilizer_to_binary_stabilizer(pauli_terms) + qvmstab = QVM_Stabilizer(num_qubits=num_qubits) + qvmstab.tableau[num_qubits:, :] = stab_mat + exp_on_i = qvmstab._rowsum_phase_accumulator(2, 3) + assert exp_on_i == 2 + + num_qubits = 2 + pauli_terms = [sZ(0) * sI(1), sZ(0) * sZ(1)] + stab_mat = pauli_stabilizer_to_binary_stabilizer(pauli_terms) + qvmstab = QVM_Stabilizer(num_qubits=num_qubits) + qvmstab.tableau[num_qubits:, :] = stab_mat + exp_on_i = qvmstab._rowsum_phase_accumulator(2, 3) + assert exp_on_i == 0 + + # now try generating random valid elements from 2-qubit group + for _ in range(100): + num_qubits = 6 + pauli_terms = [] + for _ in range(num_qubits): + pauli_terms.append(reduce(lambda x, y: x * y, [pauli_subgroup[x](idx) for idx, x in enumerate(np.random.randint(1, 4, num_qubits))])) + stab_mat = pauli_stabilizer_to_binary_stabilizer(pauli_terms) + qvmstab = QVM_Stabilizer(num_qubits=num_qubits) + qvmstab.tableau[num_qubits:, :] = stab_mat + p_on_i = qvmstab._rowsum_phase_accumulator(num_qubits, num_qubits + 1) + coeff_test = pauli_terms[1] * pauli_terms[0] + assert np.isclose(coeff_test.coefficient, (1j) ** p_on_i) + + +def test_rowsum_full(): + """ + Test the rowsum subroutine. + + This routine replaces row h with P(i) * P(h) where P is the pauli operator represented by the row + given by the index. It uses the summation to accumulate whether the phase is +1 or -1 and then + xors together the elements of the rows replacing row h + """ + # now try generating random valid elements from 2-qubit group + for _ in range(100): + num_qubits = 6 + pauli_terms = [] + for _ in range(num_qubits): + pauli_terms.append(reduce(lambda x, y: x * y, [pauli_subgroup[x](idx) for idx, x in enumerate(np.random.randint(1, 4, num_qubits))])) + stab_mat = pauli_stabilizer_to_binary_stabilizer(pauli_terms) + qvmstab = QVM_Stabilizer(num_qubits=num_qubits) + qvmstab.tableau[num_qubits:, :] = stab_mat + p_on_i = qvmstab._rowsum_phase_accumulator(num_qubits, num_qubits + 1) + if p_on_i not in [1, 3]: + coeff_test = pauli_terms[1] * pauli_terms[0] + assert np.isclose(coeff_test.coefficient, (1j) ** p_on_i) + qvmstab._rowsum(num_qubits, num_qubits + 1) + + pauli_op = binary_stabilizer_to_pauli_stabilizer(qvmstab.tableau[[num_qubits], :])[0] + true_pauli_op = pauli_terms[1] * pauli_terms[0] + assert pauli_op == true_pauli_op + + +def test_simulation_hadamard(): + """ + Test if Hadamard is applied correctly to the tableau + + The first example will be a 1 qubit example. where we perform H | 0 > = |0> + |1>. + Therefore we expect the tableau to represent the stablizer X + """ + prog = Program().inst(H(0)) + qvmstab = QVM_Stabilizer(num_qubits=1) + qvmstab._apply_hadamard(prog.instructions[0]) + x_stabilizer = np.array([[0, 1, 0], + [1, 0, 0]]) + assert np.allclose(x_stabilizer, qvmstab.tableau) + + qvmstab = QVM_Stabilizer(num_qubits=2) + qvmstab._apply_hadamard(prog.instructions[0]) + x_stabilizer = np.array([[0, 0, 1, 0, 0], + [0, 1, 0, 0, 0], + [1, 0, 0, 0, 0], + [0, 0, 0, 1, 0]]) + assert np.allclose(x_stabilizer, qvmstab.tableau) + + prog = Program().inst(H(1)) + qvmstab = QVM_Stabilizer(num_qubits=2) + qvmstab._apply_hadamard(prog.instructions[0]) + x_stabilizer = np.array([[1, 0, 0, 0, 0], + [0, 0, 0, 1, 0], + [0, 0, 1, 0, 0], + [0, 1, 0, 0, 0]]) + assert np.allclose(x_stabilizer, qvmstab.tableau) + + prog = Program().inst([H(0), H(1)]) + qvmstab = QVM_Stabilizer(num_qubits=2) + qvmstab._apply_hadamard(prog.instructions[0]) + qvmstab._apply_hadamard(prog.instructions[1]) + x_stabilizer = np.array([[0, 0, 1, 0, 0], + [0, 0, 0, 1, 0], + [1, 0, 0, 0, 0], + [0, 1, 0, 0, 0]]) + assert np.allclose(x_stabilizer, qvmstab.tableau) + + +def test_simulation_phase(): + """ + Test if the Phase gate is applied correctly to the tableau + + S|0> = |0> + """ + prog = Program().inst(S(0)) + qvmstab = QVM_Stabilizer(num_qubits=1) + qvmstab._apply_phase(prog.instructions[0]) + true_stab = np.array([[1, 1, 0], + [0, 1, 0]]) + assert np.allclose(true_stab, qvmstab.tableau) + + +def test_simulation_cnot(): + """ + Test if the simulation of CNOT is accurate + :return: + """ + prog = Program().inst([H(0), CNOT(0, 1)]) + qvmstab = QVM_Stabilizer(num_qubits=2) + qvmstab._apply_hadamard(prog.instructions[0]) + qvmstab._apply_cnot(prog.instructions[1]) + + # assert that ZZ XX stabilizes a bell state + true_stabilizers = [sX(0) * sX(1), sZ(0) * sZ(1)] + test_paulis = binary_stabilizer_to_pauli_stabilizer(qvmstab.tableau[2:, :]) + for idx, term in enumerate(test_paulis): + assert term == true_stabilizers[idx] + + +def test_measurement(): + """ + Test measurements of stabilizer state from tableau + """ + # generate bell state then measure each qubit sequentially + prog = Program().inst([H(0), CNOT(0, 1)]).measure(0, 0).measure(1, 1) + qvmstab = QVM_Stabilizer() + qvmstab.run(prog) + +if __name__ == "__main__": + test_initialization() + test_row_sum_sub_algorithm_g_test() + test_stabilizer_to_matrix_conversion() + test_rowsum_phase_accumulator() + test_rowsum_full() + test_simulation_hadamard() + test_simulation_phase() + test_simulation_cnot() + test_measurement() diff --git a/referenceqvm/unitary_generator.py b/referenceqvm/unitary_generator.py index 1cb5104..90e1733 100644 --- a/referenceqvm/unitary_generator.py +++ b/referenceqvm/unitary_generator.py @@ -418,7 +418,9 @@ def value_get(param_obj): return param_obj.index elif isinstance(param_obj, Addr): return param_obj.address - elif isinstance(param_obj, Slot): - return param_obj.value() + # elif isinstance(param_obj, Slot): + # return param_obj.value() elif isinstance(param_obj, Label): return param_obj.name + else: + raise TypeError("Object not recognizable") From 67947198d734924741afb8f8947244dd3da05a5b Mon Sep 17 00:00:00 2001 From: ncrubin Date: Fri, 27 Jul 2018 17:21:24 -0700 Subject: [PATCH 03/10] WIP stabilizer simulation Seems like we have measurement in correctly now! --- referenceqvm/qam.py | 17 +- referenceqvm/qvm_stabilizer.py | 223 ++++++++++++++++++---- referenceqvm/tests/test_stabilzier_qvm.py | 66 ++++++- 3 files changed, 256 insertions(+), 50 deletions(-) diff --git a/referenceqvm/qam.py b/referenceqvm/qam.py index c9e3507..d73b31b 100644 --- a/referenceqvm/qam.py +++ b/referenceqvm/qam.py @@ -18,6 +18,7 @@ QAM superclass. Implements state machine model, program loading and processing, kernel - leaving details of evolution up to subclasses. """ +import sys import numpy as np from referenceqvm.unitary_generator import value_get @@ -85,14 +86,15 @@ def load_program(self, pyquil_program): if not (instr.name in self.gate_set.keys() or instr.name in self.defgate_set.keys()): invalid = True break + elif isinstance(instr, Measurement): + pass else: invalid = True break # NOTE: all_inst is set by the subclass if invalid is True and self.all_inst is False: - raise TypeError("In QVM_Unitary, only Gates and DefGates are " - "supported") + raise TypeError("Some gates used are not allowed in this QAM") # set internal program and counter to their appropriate values self.program = pyquil_program @@ -105,6 +107,17 @@ def load_program(self, pyquil_program): self.num_qubits = q_max self.classical_memory = np.zeros(c_max).astype(bool) + def memory_reset(self): + if self.program is None: + raise TypeError("Program must be loaded to call reset") + + # setup quantum and classical memory + q_max, c_max = self.identify_bits() + if c_max <= 512: # allocate at least 512 cbits (as floor) + c_max = 512 + self.num_qubits = q_max + self.classical_memory = np.zeros(c_max).astype(bool) + def identify_bits(self): """ Iterates through QAM program and finds number of qubits and cbits diff --git a/referenceqvm/qvm_stabilizer.py b/referenceqvm/qvm_stabilizer.py index 5e6c39e..0ca33a8 100644 --- a/referenceqvm/qvm_stabilizer.py +++ b/referenceqvm/qvm_stabilizer.py @@ -20,7 +20,7 @@ """ import sys from functools import reduce -from pyquil.quil import Program +from pyquil.quil import Program, get_classical_addresses_from_program from pyquil.quilbase import * from pyquil.paulis import PauliTerm, sI, sZ, sX, sY @@ -56,12 +56,69 @@ def __init__(self, num_qubits=None, program=None, program_counter=None, classical_memory=classical_memory, gate_set=gate_set, defgate_set=defgate_set) - + # this sets that arbitrary instructions are not allowed + # this can probably be factored out eventually. It is a relic of the + # the old days of generating unitaries of programs + self.all_inst = False if num_qubits is None: self.tableau = None else: self.tableau = self._n_qubit_tableau(num_qubits) + def load_program(self, pyquil_program): + """ + Loads a pyQuil program into the QAM memory. + Synthesizes pyQuil program into instruction list. + Initializes program object and program counter. + + This overwrites the parent class load program. Required because + unitary check QVM checks makes it annoying to allow measurements. + TODO: Change load_program for each qvm-subclass + + :param Program pyquil_program: program to be ran + """ + # typecheck + if not isinstance(pyquil_program, Program): + raise TypeError("I can only generate from pyQuil programs") + + if self.all_inst is None: + raise NotImplementedError("QAM needs to be subclassed in order to " + "load program") + + # create defgate dictionary + defined_gates = {} + for dg in pyquil_program.defined_gates: + defined_gates[dg.name] = dg.matrix + self.defgate_set = defined_gates + + # if QVM_Unitary, check if all instructions are valid. + invalid = False + for instr in pyquil_program: + if isinstance(instr, Gate): + if not (instr.name in self.gate_set.keys() or instr.name in self.defgate_set.keys()): + invalid = True + break + elif isinstance(instr, Measurement): + pass + else: + invalid = True + break + + # NOTE: all_inst is set by the subclass + if invalid is True and self.all_inst is False: + raise TypeError("Some gates used are not allowed in this QAM") + + # set internal program and counter to their appropriate values + self.program = pyquil_program + self.program_counter = 0 + + # setup quantum and classical memory + q_max, c_max = self.identify_bits() + if c_max <= 512: # allocate at least 512 cbits (as floor) + c_max = 512 + self.num_qubits = q_max + self.classical_memory = np.zeros(c_max).astype(bool) + def _n_qubit_tableau(self, num_qubits): """ Construct an empty tableau for a n-qubit system @@ -77,17 +134,98 @@ def _n_qubit_tableau(self, num_qubits): tableau[ii, ii] = 1 return tableau - def run(self, pyquil_program): + def transition(self, instruction): + """ + Implements a full transition, including pre/post noise hooks. + + :param QuilAction instruction: instruction to be executed + + :return: if program is halted TRUE, else FALSE + :rtype: bool int + """ + self.pre() + self._transition(instruction) + self.post() + + # return HALTED (i.e. program_counter is end of program) + return self.program_counter == len(self.program) + + def pre(self): + """ + Pre-transition hook - use for noisy evolution models. Unimplemented for now. + """ + pass + + def post(self): + """ + Post-transition hook - use for noisy evolution models. Unimplemented for now. + """ + pass + + def _transition(self, instruction): """ - Run comman + Implements a transition on the generator matrix representing the stabilizers + + :param Gate instruction: QuilAction gate to be implemented + """ + if isinstance(instruction, Measurement): + # mutates classical memory! I should change this... + self._apply_measurement(instruction) + self.program_counter += 1 + + elif isinstance(instruction, Gate): + # apply Gate or DefGate + if instruction.name == 'H': + self._apply_hadamard(instruction) + elif instruction.name == 'S': + self._apply_phase(instruction) + elif instruction.name == 'CNOT': + self._apply_cnot(instruction) + elif instruction.name == 'I': + pass + else: + raise TypeError("We checked for correct gate types previously" + + " so the impossible has happened!") + + self.program_counter += 1 + + elif isinstance(instruction, Jump): + # unconditional Jump; go directly to Label + self.program_counter = self.find_label(instruction.target) + + elif isinstance(instruction, JumpTarget): + # Label; pass straight over + self.program_counter += 1 + + def run(self, pyquil_program, classical_addresses=None, trials=1): + """ + Run program. + + Loads and checks program if all gates are within the stabilizer set. + Then executes program :param porgram: + :param classical_addresses: + :param trials: :return: """ self.load_program(pyquil_program) - # set up stabilizers - self.tableau = self._n_qubit_tableau(self.num_qubits) + if classical_addresses is None: + classical_addresses = get_classical_addresses_from_program(pyquil_program) + + results = [] + for _ in range(trials): + # set up stabilizers + self.tableau = self._n_qubit_tableau(self.num_qubits) + self.kernel() + results.append(list(map(int, self.classical_memory[classical_addresses]))) + + # reset qvm + self.memory_reset() + self.program_counter = 0 + + return results def stabilizer_tableau(self): """ @@ -145,7 +283,6 @@ def _rowsum_phase_accumulator(self, h, i): self.tableau[i, self.num_qubits + j], self.tableau[h, j], self.tableau[h, self.num_qubits + j]) - phase_accumulator += 2 * self.tableau[h, -1] phase_accumulator += 2 * self.tableau[i, -1] return phase_accumulator % 4 @@ -244,8 +381,9 @@ def _apply_measurement(self, instruction): t_qbit = value_get(instruction.qubit) t_cbit = value_get(instruction.classical_reg) - print(self.tableau) - check_xpa = False # check if x_pa = 1 for p in {n + 1 ... 2 n} of tableau + # check if the output of the measurement is random + # this is analogous to the case when the measurement operator does not + # commute with at least one stabilizer if any(self.tableau[self.num_qubits:, t_qbit] == 1): # find the first one. xpa_idx = np.where(self.tableau[self.num_qubits:, t_qbit] == 1)[0][0] # take the first index @@ -253,21 +391,36 @@ def _apply_measurement(self, instruction): for ii in range(2 * self.num_qubits): # loop over each row and call rowsum(ii, xpa_idx) if ii != xpa_idx and self.tableau[ii, t_qbit] == 1: self._rowsum(ii, xpa_idx) + + # moving the operator into the destabilizer and then replacing with + # the measurement operator self.tableau[xpa_idx - self.num_qubits, :] = self.tableau[xpa_idx, :] + # this is like replacing the non-commuting element with the measurement operator self.tableau[xpa_idx, :] = np.zeros((1, 2 * self.num_qubits + 1)) self.tableau[xpa_idx, t_qbit + self.num_qubits] = 1 + # perform the measurement self.tableau[xpa_idx, -1] = 1 if np.random.random() > 0.5 else 0 # set classical results to return self.classical_memory[t_cbit] = self.tableau[xpa_idx, -1] + # outcome of measurement is deterministic...need to determine sign else: # augment tableaue with a scratch space - self.tableau = np.vstack((self.tableau, np.zeros((1, 2 * self.num_qubits + 1)))) + self.tableau = np.vstack((self.tableau, np.zeros((1, 2 * self.num_qubits + 1), dtype=int))) for ii in range(self.num_qubits): - self._rowsum(2 * self.num_qubits + 1, ii) + # We need to check if R(i) anticommutes with Za...which it does if x_{ia} = 1 + if self.tableau[ii, t_qbit] == 1: # refrencing the destabilizers + + # check something silly. Does the destabilizer anticommute with the observable? It SHOULD! + tmp_vector_representing_z_qubit = np.zeros((2 * self.num_qubits), dtype=int) + tmp_vector_representing_z_qubit[t_qbit] = 1 + assert symplectic_inner_product(tmp_vector_representing_z_qubit, self.tableau[ii, :-1]) == 1 + + # row sum on the stabilizers (summing up operators such that we get operator Z_{a}) + self._rowsum(2 * self.num_qubits, ii + self.num_qubits) # note: A-G says 2 n + 1...this is correct...but they start counting at 1 not zero # set the classical bit to be the last element of the scratch row self.classical_memory[t_cbit] = self.tableau[-1, -1] @@ -275,35 +428,6 @@ def _apply_measurement(self, instruction): # remove the scratch space self.tableau = self.tableau[:2 * self.num_qubits, :] - def transition(self, instruction): - """ - Implements a transition on the generator matrix representing the stabilizers - - :param Gate instruction: QuilAction gate to be implemented - """ - if isinstance(instruction, Measurement): - # mutates classical memory! I should change this... - self._apply_measurement(instruction) - self.program_counter += 1 - - elif isinstance(instruction, Gate): - # apply Gate or DefGate - if Gate.name == 'H': - self._apply_hadamard(instruction) - elif Gate.name == 'S': - self._apply_phase(instruction) - elif Gate.name == 'CNOT': - self._apply_cnot(instruction) - self.program_counter += 1 - - elif isinstance(instruction, Jump): - # unconditional Jump; go directly to Label - self.program_counter = self.find_label(instruction.target) - - elif isinstance(instruction, JumpTarget): - # Label; pass straight over - self.program_counter += 1 - def pauli_stabilizer_to_binary_stabilizer(stabilizer_list): """ @@ -414,6 +538,25 @@ def binary_rref(code_matrix): return code_matrix +def symplectic_inner_product(vector1, vector2): + """ + Operators commute if their binary form symplectic inner product is zero + + Operators anticommute if their binary form symplectic inner product is one + + :param vector1: binary form of operator with no sign info + :param vector2: binary form of a pauli operator with no sign info + :return: 0, 1 + """ + if vector1.shape != vector2.shape: + raise ValueError("vectors must be the same size.") + + # TODO: add a check for binary or integer linear arrays + + hadamard_product = np.multiply(vector1, vector2) + return reduce(lambda x, y: x ^ y, hadamard_product) + + if __name__ == "__main__": # practice my clifford transformations # = |0000...0N> diff --git a/referenceqvm/tests/test_stabilzier_qvm.py b/referenceqvm/tests/test_stabilzier_qvm.py index 0a79728..9322abf 100644 --- a/referenceqvm/tests/test_stabilzier_qvm.py +++ b/referenceqvm/tests/test_stabilzier_qvm.py @@ -7,7 +7,7 @@ from functools import reduce from pyquil.paulis import sI, sX, sY, sZ, PAULI_COEFF, PAULI_OPS from pyquil.quil import Program -from pyquil.gates import H, S, CNOT +from pyquil.gates import H, S, CNOT, I from referenceqvm.qvm_stabilizer import (QVM_Stabilizer, pauli_stabilizer_to_binary_stabilizer, binary_stabilizer_to_pauli_stabilizer) from referenceqvm.unitary_generator import value_get @@ -183,7 +183,12 @@ def test_rowsum_full(): pauli_terms = [] for _ in range(num_qubits): pauli_terms.append(reduce(lambda x, y: x * y, [pauli_subgroup[x](idx) for idx, x in enumerate(np.random.randint(1, 4, num_qubits))])) - stab_mat = pauli_stabilizer_to_binary_stabilizer(pauli_terms) + try: + stab_mat = pauli_stabilizer_to_binary_stabilizer(pauli_terms) + except: + # we have to do this because I'm not strictly making valid n-qubit + # stabilizers + continue qvmstab = QVM_Stabilizer(num_qubits=num_qubits) qvmstab.tableau[num_qubits:, :] = stab_mat p_on_i = qvmstab._rowsum_phase_accumulator(num_qubits, num_qubits + 1) @@ -191,8 +196,10 @@ def test_rowsum_full(): coeff_test = pauli_terms[1] * pauli_terms[0] assert np.isclose(coeff_test.coefficient, (1j) ** p_on_i) qvmstab._rowsum(num_qubits, num_qubits + 1) - - pauli_op = binary_stabilizer_to_pauli_stabilizer(qvmstab.tableau[[num_qubits], :])[0] + try: + pauli_op = binary_stabilizer_to_pauli_stabilizer(qvmstab.tableau[[num_qubits], :])[0] + except: + continue true_pauli_op = pauli_terms[1] * pauli_terms[0] assert pauli_op == true_pauli_op @@ -270,14 +277,54 @@ def test_simulation_cnot(): assert term == true_stabilizers[idx] -def test_measurement(): +def test_measurement_noncommuting(): """ Test measurements of stabilizer state from tableau + + This tests the non-commuting measurement operator """ # generate bell state then measure each qubit sequentially - prog = Program().inst([H(0), CNOT(0, 1)]).measure(0, 0).measure(1, 1) + prog = Program().inst(H(0)).measure(0, 0) qvmstab = QVM_Stabilizer() - qvmstab.run(prog) + results = qvmstab.run(prog, trials=5000) + assert np.isclose(np.mean(results), 0.5, rtol=0.1) + + +def test_measurement_commuting(): + """ + Test measuremt of stabilzier state from tableau + + This tests when the measurement operator commutes with the stabilizers + + To test we will first test draw's from a blank stabilizer and then with X + """ + identity_program = Program().inst([I(0)]).measure(0, 0) + qvmstab = QVM_Stabilizer() + results = qvmstab.run(identity_program, trials=1000) + assert all(np.array(results) == 0) + + +def test_measurement_commuting_result_one(): + """ + Test measuremt of stabilzier state from tableau + + This tests when the measurement operator commutes with the stabilizers + + This time we will generate the stabilizer -Z|1> = |1> so we we need to do + a bitflip...not just identity. A Bitflip is HSSH = X + """ + identity_program = Program().inst([H(0), S(0), S(0), H(0)]).measure(0, 0) + qvmstab = QVM_Stabilizer() + results = qvmstab.run(identity_program, trials=1000) + assert all(np.array(results) == 1) + + +def test_bell_state_measurements(): + prog = Program().inst(H(0), CNOT(0, 1)).measure(0, 0).measure(1, 1) + qvmstab = QVM_Stabilizer() + results = qvmstab.run(prog, trials=5000) + assert np.isclose(np.mean(results), 0.5, rtol=0.1) + if __name__ == "__main__": test_initialization() @@ -288,4 +335,7 @@ def test_measurement(): test_simulation_hadamard() test_simulation_phase() test_simulation_cnot() - test_measurement() + test_measurement_noncommuting() + test_measurement_commuting() + test_measurement_commuting_result_one() + test_bell_state_measurements() From e97e812bf3802c5e4bccb10875e8c733ff42307c Mon Sep 17 00:00:00 2001 From: ncrubin Date: Mon, 30 Jul 2018 12:58:33 -0700 Subject: [PATCH 04/10] Stabilizer sim can now project out state and rho Stabilizer simulator can now return the state (which has been tested) and the density matrix (which has not been tested) More tests need to be generated to validate that the stabilizer simulation actually works. --- referenceqvm/qvm_stabilizer.py | 25 ++- referenceqvm/state_actions.py | 87 +++++++- referenceqvm/tests/test_state_actions.py | 268 +++++++++++++++++++++++ 3 files changed, 377 insertions(+), 3 deletions(-) create mode 100644 referenceqvm/tests/test_state_actions.py diff --git a/referenceqvm/qvm_stabilizer.py b/referenceqvm/qvm_stabilizer.py index 0ca33a8..9e710ca 100644 --- a/referenceqvm/qvm_stabilizer.py +++ b/referenceqvm/qvm_stabilizer.py @@ -16,7 +16,8 @@ ############################################################################## """ Pure QVM that only executes pyQuil programs containing Clifford group generators, -and return the wavefunction or stabilizer +and return the wavefunction or stabilizer. This is based off of the paper by +Aaronson and Gottesman Phys. Rev. A 70, 052328 """ import sys from functools import reduce @@ -26,7 +27,7 @@ from referenceqvm.qam import QAM from referenceqvm.gates import stabilizer_gate_matrix -from referenceqvm.unitary_generator import value_get +from referenceqvm.unitary_generator import value_get, tensor_up class QVM_Stabilizer(QAM): @@ -227,6 +228,25 @@ def run(self, pyquil_program, classical_addresses=None, trials=1): return results + def density(self, pyquil_program): + """ + Run program and compute the density matrix + + Loads and checks program if all gates are within the stabilizer set. + Then executes program + + :param porgram: + :return: + """ + self.load_program(pyquil_program) + self.tableau = self._n_qubit_tableau(self.num_qubits) + self.kernel() + stabilizers = binary_stabilizer_to_pauli_stabilizer(self.stabilizer_tableau()) + pauli_ops = reduce(lambda x, y: x * (1 + y), stabilizers) + pauli_ops *= (2 ** (-self.num_qubits)) + + return tensor_up(pauli_ops, self.num_qubits) + def stabilizer_tableau(self): """ return the stabilizer part of the tableau @@ -429,6 +449,7 @@ def _apply_measurement(self, instruction): self.tableau = self.tableau[:2 * self.num_qubits, :] + def pauli_stabilizer_to_binary_stabilizer(stabilizer_list): """ Convert a list of stabilizers represented as PauliTerms to a binary tableau form diff --git a/referenceqvm/state_actions.py b/referenceqvm/state_actions.py index 2e85709..a036f80 100644 --- a/referenceqvm/state_actions.py +++ b/referenceqvm/state_actions.py @@ -9,6 +9,7 @@ Given """ +from scipy.sparse import csc_matrix import numpy as np from pyquil.paulis import PauliTerm @@ -17,10 +18,14 @@ def compute_action(classical_state, pauli_operator, num_qubits): """ Compute action of Pauli opertors on a classical state + The classical state is enumerated as the left most bit is the least-significant + bit. This is how one usually reads off classical data from the QVM. Not + how the QVM stores computational basis states. + :param classical_state: binary repr of a state or an integer. Should be left most bit (0th position) is the most significant bit :param num_qubits: - :return: + :return: new classical state and the coefficient it picked up. """ if not isinstance(pauli_operator, PauliTerm): raise TypeError("pauli_operator must be a PauliTerm") @@ -37,3 +42,83 @@ def compute_action(classical_state, pauli_operator, num_qubits): if len(classical_state) != num_qubits: raise TypeError("classical state not long enough") + # iterate through tensor elements of pauli_operator + new_classical_state = classical_state.copy() + coefficient = 1 + for qidx, telem in pauli_operator: + if telem == 'X': + new_classical_state[qidx] = new_classical_state[qidx] ^ 1 + elif telem == 'Y': + new_classical_state[qidx] = new_classical_state[qidx] ^ 1 + # set coeff + if new_classical_state[qidx] == 0: + coefficient *= -1j + else: + coefficient *= 1j + + elif telem == 'Z': + # set coeff + if new_classical_state[qidx] == 1: + coefficient *= -1 + + return new_classical_state, coefficient + + +def state_family_generator(state, pauli_operator): + """ + Generate a new state by applying the pauli_operator to each computational bit-string + + This is accomplished in a sparse format where a sparse vector is returned + after the action is accumulate in a new list of data and indices + + :param state: + :param pauli_operator: + :return: + """ + if not isinstance(state, csc_matrix): + raise TypeError("we only take csc_matrix") + + num_qubits = int(np.log2(state.shape[0])) + new_coeffs = [] + new_indices = [] + + # iterate through non-zero + rindices, cindices = state.nonzero() + for ridx, cidx in zip(rindices, cindices): + # this is so gross looking + bitstring = list(map(int, np.binary_repr(ridx, width=num_qubits)))[::-1] + new_ket, new_coefficient = compute_action(bitstring, pauli_operator, num_qubits) + new_indices.append(int("".join([str(x) for x in new_ket[::-1]]), 2)) + new_coeffs.append(state[ridx, cidx] * new_coefficient) + + return csc_matrix((new_coeffs, (new_indices, [0] * len(new_indices))), + shape=(2 ** num_qubits, 1), dtype=complex) + + +def project_stabilized_state(stabilizer_list, num_qubits=None): + """ + Project out the state stabilized by the stabilizer matrix + + |psi> = (1/2^{n}) * Product_{i=0}{n-1}[ 1 + G_{i}] |vac> + + :param stabilizer_list: + :return: + """ + if num_qubits is None: + num_qubits = len(stabilizer_list) + + indptr = np.array([0]) + indices = np.array([0]) + data = np.array([1]) + state = csc_matrix((data, (indices, indptr)), shape=(2 ** num_qubits, 1), + dtype=complex) + + for generator in stabilizer_list: + # (I + G(i)) / 2 + state += state_family_generator(state, generator) + state /= 2 + + normalization = (state.conj().T.dot(state)).todense() + state /= np.sqrt(float(normalization.real)) # this is needed or it will cast as a np.matrix + + return state diff --git a/referenceqvm/tests/test_state_actions.py b/referenceqvm/tests/test_state_actions.py new file mode 100644 index 0000000..9240512 --- /dev/null +++ b/referenceqvm/tests/test_state_actions.py @@ -0,0 +1,268 @@ +""" +Test the infrastructure for building a state by projection onto the +1 +eigenspace of a set of generators or stabilizers +""" +import numpy as np +from referenceqvm.state_actions import compute_action, project_stabilized_state +from referenceqvm.tests.test_stabilzier_qvm import five_qubit_code_generators +from pyquil.paulis import sX, sZ, sY, sI, PauliSum +import pytest + + +def test_compute_action_type_checks(): + """ + Make sure type checks are consistent and working + """ + with pytest.raises(TypeError): + compute_action([0, 0, 0, 0, 0], PauliSum([sX(0)]), 5) + + with pytest.raises(TypeError): + compute_action([0, 0, 0, 0, 0], sX(0), 4) + + with pytest.raises(TypeError): + compute_action(3, 'a', 4) + + with pytest.raises(TypeError): + compute_action(-3, sX(0), 4) + + with pytest.raises(TypeError): + compute_action('0001', sX(0), 4) + + +def test_compute_action_identity(): + """ + Action of Pauli operators on state + """ + comp_basis_state = [0, 0, 0, 0] + for ii in range(4): + pauli_term = sI(ii) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + assert new_basis_state == comp_basis_state + assert np.isclose(coeff, 1) + + +def test_compute_action_X(): + """ + Action of Pauli operators on state + """ + comp_basis_state = [0, 0, 0, 0] + for ii in range(4): + pauli_term = sX(ii) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] = 1 + assert new_basis_state == true_basis_state + assert np.isclose(coeff, 1) + + comp_basis_state = [1, 1, 1, 1] + for ii in range(4): + pauli_term = sX(ii) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] = true_basis_state[ii] ^ 1 + assert new_basis_state == true_basis_state + assert np.isclose(coeff, 1) + + +def test_compute_action_XX(): + """ + Action of Pauli operators on state + """ + comp_basis_state = [0, 0, 0, 0] + for ii in range(3): + pauli_term = sX(ii) * sX(ii + 1) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] = true_basis_state[ii + 1] = 1 + assert new_basis_state == true_basis_state + assert np.isclose(coeff, 1) + + comp_basis_state = [1, 1, 1, 1] + for ii in range(3): + pauli_term = sX(ii) * sX(ii + 1) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] = true_basis_state[ii] ^ 1 + true_basis_state[ii + 1] = true_basis_state[ii + 1] ^ 1 + assert new_basis_state == true_basis_state + assert np.isclose(coeff, 1) + + +def test_compute_action_Y(): + """ + Action of Pauli operators on state + """ + comp_basis_state = [0, 0, 0, 0] + for ii in range(4): + pauli_term = sY(ii) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] = true_basis_state[ii] ^ 1 + assert new_basis_state == true_basis_state + assert np.isclose(coeff, 1j) + + comp_basis_state = [1, 1, 1, 1] + for ii in range(4): + pauli_term = sY(ii) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] = true_basis_state[ii] ^ 1 + assert new_basis_state == true_basis_state + assert np.isclose(coeff, -1j) + + +def test_compute_action_YY(): + """ + Action of Pauli operators on state + """ + comp_basis_state = [0, 0, 0, 0] + for ii in range(3): + pauli_term = sY(ii) * sY(ii + 1) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] = true_basis_state[ii + 1] = 1 + assert new_basis_state == true_basis_state + assert np.isclose(coeff, -1) + + comp_basis_state = [1, 1, 1, 1] + for ii in range(3): + pauli_term = sY(ii) * sY(ii + 1) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] = true_basis_state[ii] ^ 1 + true_basis_state[ii + 1] = true_basis_state[ii + 1] ^ 1 + assert new_basis_state == true_basis_state + assert np.isclose(coeff, -1) + + +def test_compute_action_Z(): + """ + Action of Pauli operators on state + """ + comp_basis_state = [0, 0, 0, 0] + for ii in range(4): + pauli_term = sZ(ii) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] = true_basis_state[ii] + assert new_basis_state == true_basis_state + assert np.isclose(coeff, 1) + + comp_basis_state = [1, 1, 1, 1] + for ii in range(4): + pauli_term = sZ(ii) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] = true_basis_state[ii] + assert new_basis_state == true_basis_state + assert np.isclose(coeff, -1) + + +def test_compute_action_ZZ(): + """ + Action of Pauli operators on state + """ + comp_basis_state = [0, 0, 0, 0] + for ii in range(3): + pauli_term = sZ(ii) * sZ(ii + 1) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + assert new_basis_state == true_basis_state + assert np.isclose(coeff, 1) + + comp_basis_state = [1, 1, 1, 1] + for ii in range(3): + pauli_term = sZ(ii) * sZ(ii + 1) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + assert new_basis_state == true_basis_state + assert np.isclose(coeff, 1) + + +def test_compute_action_XY(): + """ + Action of Pauli operators on state + """ + comp_basis_state = [0, 0, 0, 0] + for ii in range(3): + pauli_term = sX(ii) * sY(ii + 1) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] ^= 1 + true_basis_state[ii + 1] ^= 1 + assert new_basis_state == true_basis_state + assert np.isclose(coeff, 1j) + + comp_basis_state = [1, 1, 1, 1] + for ii in range(3): + pauli_term = sX(ii) * sY(ii + 1) + new_basis_state, coeff = compute_action(comp_basis_state, pauli_term, + len(comp_basis_state)) + # abuse of comparisons in python + true_basis_state = comp_basis_state.copy() + true_basis_state[ii] ^= 1 + true_basis_state[ii + 1] ^= 1 + assert new_basis_state == true_basis_state + assert np.isclose(coeff, -1j) + + +def test_stabilizer_projection_Z(): + """ + test if we project out the correct state + """ + stabilizer_state = project_stabilized_state([sZ(0)]) + true_state = np.zeros((2, 1)) + true_state[0, 0] = 1 + assert np.allclose(true_state, stabilizer_state.todense()) + + +def test_stabilizer_projection_ZZ(): + """ + test if we project out the correct state + """ + stabilizer_state = project_stabilized_state([sZ(0) * sZ(1), sX(0) * sX(1)]) + true_state = np.zeros((4, 1)) + true_state[0, 0] = true_state[3, 0] = 1 + true_state /= np.sqrt(2) + assert np.allclose(true_state, stabilizer_state.todense()) + + +def test_stabilizer_projection_ZZZ(): + """ + test if we project out the correct state + """ + stabilizer_state = project_stabilized_state([sZ(0) * sZ(1), sZ(1) * sZ(2), + sX(0) * sX(1) * sX(2)]) + true_state = np.zeros((8, 1)) + true_state[0, 0] = true_state[7, 0] = 1 + true_state /= np.sqrt(2) + assert np.allclose(true_state, stabilizer_state.todense()) + From 4a07a099755290d7375eb890997d9904dc60232b Mon Sep 17 00:00:00 2001 From: ncrubin Date: Tue, 31 Jul 2018 11:29:14 -0700 Subject: [PATCH 05/10] Working tested Improved Stabilizer Simulation Error was on how we were pulling the qubits from the instruction. NOTE: Do not use get_qubits() as it always return qubit indices in numerical order. --- referenceqvm/qam.py | 2 - referenceqvm/qvm_stabilizer.py | 31 ++++-- referenceqvm/state_actions.py | 30 ++++-- referenceqvm/tests/test_stabilzier_qvm.py | 124 +++++++++++++++++++--- referenceqvm/tests/test_state_actions.py | 4 +- 5 files changed, 162 insertions(+), 29 deletions(-) diff --git a/referenceqvm/qam.py b/referenceqvm/qam.py index d73b31b..05a71ee 100644 --- a/referenceqvm/qam.py +++ b/referenceqvm/qam.py @@ -86,8 +86,6 @@ def load_program(self, pyquil_program): if not (instr.name in self.gate_set.keys() or instr.name in self.defgate_set.keys()): invalid = True break - elif isinstance(instr, Measurement): - pass else: invalid = True break diff --git a/referenceqvm/qvm_stabilizer.py b/referenceqvm/qvm_stabilizer.py index 9e710ca..00bbd1b 100644 --- a/referenceqvm/qvm_stabilizer.py +++ b/referenceqvm/qvm_stabilizer.py @@ -20,14 +20,17 @@ Aaronson and Gottesman Phys. Rev. A 70, 052328 """ import sys +import numpy as np from functools import reduce from pyquil.quil import Program, get_classical_addresses_from_program from pyquil.quilbase import * from pyquil.paulis import PauliTerm, sI, sZ, sX, sY +from pyquil.wavefunction import Wavefunction from referenceqvm.qam import QAM from referenceqvm.gates import stabilizer_gate_matrix from referenceqvm.unitary_generator import value_get, tensor_up +from referenceqvm.state_actions import project_stabilized_state class QVM_Stabilizer(QAM): @@ -242,11 +245,27 @@ def density(self, pyquil_program): self.tableau = self._n_qubit_tableau(self.num_qubits) self.kernel() stabilizers = binary_stabilizer_to_pauli_stabilizer(self.stabilizer_tableau()) - pauli_ops = reduce(lambda x, y: x * (1 + y), stabilizers) - pauli_ops *= (2 ** (-self.num_qubits)) - + pauli_ops = list(map(lambda x: 0.5 * (sI(0) + x), stabilizers)) + pauli_ops = reduce(lambda x, y: x * y, pauli_ops) return tensor_up(pauli_ops, self.num_qubits) + def wavefunction(self, program): + """ + Simulate the stabilizer program and then project out the final state. + + Return the final state as a wavefunction pyquil object + + :param program: pyQuil program composed of stabilizer operations only + :return: pyquil.Wavefunction.wavefunction object. + """ + self.load_program(program) + self.tableau = self._n_qubit_tableau(self.num_qubits) + self.kernel() + stabilizers = binary_stabilizer_to_pauli_stabilizer(self.stabilizer_tableau()) + stabilizer_state = project_stabilized_state(stabilizers) + stabilizer_state = np.array(stabilizer_state.todense()) # todense() returns a matrixlib type + return Wavefunction(stabilizer_state.flatten()) + def stabilizer_tableau(self): """ return the stabilizer part of the tableau @@ -355,7 +374,7 @@ def _apply_cnot(self, instruction): :param instruction: pyquil abstract instruction. Must have """ - a, b = list(instruction.get_qubits()) # control (a) and target (b) + a, b = [value_get(x) for x in instruction.qubits] for i in range(2 * self.num_qubits): self.tableau[i, -1] = self._cnot_phase_update(i, a, b) self.tableau[i, b] = self.tableau[i, b] ^ self.tableau[i, a] @@ -380,7 +399,7 @@ def _apply_hadamard(self, instruction): :param instruction: :return: """ - qubit_label = list(instruction.get_qubits())[0] + qubit_label = [value_get(x) for x in instruction.qubits][0] for i in range(2 * self.num_qubits): self.tableau[i, -1] = self.tableau[i, -1] ^ (self.tableau[i, qubit_label] * self.tableau[i, qubit_label + self.num_qubits]) self.tableau[i, [qubit_label, qubit_label + self.num_qubits]] = self.tableau[i, [qubit_label + self.num_qubits, qubit_label]] @@ -392,7 +411,7 @@ def _apply_phase(self, instruction): :param instruction: :return: """ - qubit_label = list(instruction.get_qubits())[0] + qubit_label = [value_get(x) for x in instruction.qubits][0] for i in range(2 * self.num_qubits): self.tableau[i, -1] = self.tableau[i, -1] ^ (self.tableau[i, qubit_label] * self.tableau[i, qubit_label + self.num_qubits]) self.tableau[i, qubit_label + self.num_qubits] = self.tableau[i, qubit_label + self.num_qubits] ^ self.tableau[i, qubit_label] diff --git a/referenceqvm/state_actions.py b/referenceqvm/state_actions.py index a036f80..4253875 100644 --- a/referenceqvm/state_actions.py +++ b/referenceqvm/state_actions.py @@ -9,6 +9,7 @@ Given """ +import sys from scipy.sparse import csc_matrix import numpy as np from pyquil.paulis import PauliTerm @@ -89,27 +90,44 @@ def state_family_generator(state, pauli_operator): bitstring = list(map(int, np.binary_repr(ridx, width=num_qubits)))[::-1] new_ket, new_coefficient = compute_action(bitstring, pauli_operator, num_qubits) new_indices.append(int("".join([str(x) for x in new_ket[::-1]]), 2)) - new_coeffs.append(state[ridx, cidx] * new_coefficient) + new_coeffs.append(state[ridx, cidx] * new_coefficient * pauli_operator.coefficient) return csc_matrix((new_coeffs, (new_indices, [0] * len(new_indices))), shape=(2 ** num_qubits, 1), dtype=complex) -def project_stabilized_state(stabilizer_list, num_qubits=None): +def project_stabilized_state(stabilizer_list, num_qubits=None, + classical_state=None): """ Project out the state stabilized by the stabilizer matrix |psi> = (1/2^{n}) * Product_{i=0}{n-1}[ 1 + G_{i}] |vac> :param stabilizer_list: - :return: + :param num_qubits: integer number of qubits + :param classical_state: Default None. Defaults to |+>^{otimes n} + + :return: state projected by stabilizers """ if num_qubits is None: num_qubits = len(stabilizer_list) - indptr = np.array([0]) - indices = np.array([0]) - data = np.array([1]) + if classical_state is None: + indptr = np.array([0] * (2 ** num_qubits)) + indices = np.arange(2 ** num_qubits) + data = np.ones((2 ** num_qubits)) / np.sqrt((2 ** num_qubits)) + else: + if not isinstance(classical_state, list): + raise TypeError("I only accept lists as the classical state") + if len(classical_state) != num_qubits: + raise TypeError("Classical state does not match the number of qubits") + + # convert into an integer + ket_idx = int("".join([str(x) for x in classical_state[::-1]]), 2) + indptr = np.array([0]) + indices = np.array([ket_idx]) + data = np.array([1.]) + state = csc_matrix((data, (indices, indptr)), shape=(2 ** num_qubits, 1), dtype=complex) diff --git a/referenceqvm/tests/test_stabilzier_qvm.py b/referenceqvm/tests/test_stabilzier_qvm.py index 9322abf..c915315 100644 --- a/referenceqvm/tests/test_stabilzier_qvm.py +++ b/referenceqvm/tests/test_stabilzier_qvm.py @@ -10,7 +10,9 @@ from pyquil.gates import H, S, CNOT, I from referenceqvm.qvm_stabilizer import (QVM_Stabilizer, pauli_stabilizer_to_binary_stabilizer, binary_stabilizer_to_pauli_stabilizer) -from referenceqvm.unitary_generator import value_get +from referenceqvm.state_actions import project_stabilized_state +from referenceqvm.api import QVMConnection + pauli_subgroup = [sI, sX, sY, sZ] five_qubit_code_generators = [sX(0) * sZ(1) * sZ(2) * sX(3) * sI(4), @@ -245,12 +247,26 @@ def test_simulation_hadamard(): [0, 1, 0, 0, 0]]) assert np.allclose(x_stabilizer, qvmstab.tableau) + prog = Program().inst([H(0), H(1), H(0), H(1)]) + qvmstab = QVM_Stabilizer(num_qubits=2) + qvmstab._apply_hadamard(prog.instructions[0]) + qvmstab._apply_hadamard(prog.instructions[1]) + qvmstab._apply_hadamard(prog.instructions[2]) + qvmstab._apply_hadamard(prog.instructions[3]) + x_stabilizer = np.array([[1, 0, 0, 0, 0], + [0, 1, 0, 0, 0], + [0, 0, 1, 0, 0], + [0, 0, 0, 1, 0]]) + assert np.allclose(x_stabilizer, qvmstab.tableau) + def test_simulation_phase(): """ Test if the Phase gate is applied correctly to the tableau S|0> = |0> + + S|+> = |R> """ prog = Program().inst(S(0)) qvmstab = QVM_Stabilizer(num_qubits=1) @@ -259,6 +275,14 @@ def test_simulation_phase(): [0, 1, 0]]) assert np.allclose(true_stab, qvmstab.tableau) + prog = Program().inst([H(0), S(0)]) + qvmstab = QVM_Stabilizer(num_qubits=1) + qvmstab._apply_hadamard(prog.instructions[0]) + qvmstab._apply_phase(prog.instructions[1]) + true_stab = np.array([[0, 1, 0], + [1, 1, 0]]) + assert np.allclose(true_stab, qvmstab.tableau) + def test_simulation_cnot(): """ @@ -276,6 +300,41 @@ def test_simulation_cnot(): for idx, term in enumerate(test_paulis): assert term == true_stabilizers[idx] + # test that CNOT does nothing to |00> state + prog = Program().inst([CNOT(0, 1)]) + qvmstab = QVM_Stabilizer(num_qubits=2) + qvmstab._apply_cnot(prog.instructions[0]) + true_tableau = np.array([[1, 1, 0, 0, 0], # X1 -> X1 X2 + [0, 1, 0, 0, 0], # X2 -> X2 + [0, 0, 1, 0, 0], # Z1 -> Z1 + [0, 0, 1, 1, 0]]) # Z2 -> Z1 Z2 + + # note that Z1, Z1 Z2 still stabilizees |00> + assert np.allclose(true_tableau, qvmstab.tableau) + + + # test that CNOT produces 11 state after X + prog = Program().inst([H(0), S(0), S(0), H(0), CNOT(0, 1)]) + qvmstab = QVM_Stabilizer(num_qubits=2) + qvmstab._apply_hadamard(prog.instructions[0]) + qvmstab._apply_phase(prog.instructions[1]) + qvmstab._apply_phase(prog.instructions[2]) + qvmstab._apply_hadamard(prog.instructions[3]) + qvmstab._apply_cnot(prog.instructions[4]) + true_tableau = np.array([[1, 1, 0, 0, 0], + [0, 1, 0, 0, 0], + [0, 0, 1, 0, 1], + [0, 0, 1, 1, 0]]) + + # note that -Z1, Z1 Z2 still stabilizees |11> + assert np.allclose(true_tableau, qvmstab.tableau) + + test_paulis = binary_stabilizer_to_pauli_stabilizer(qvmstab.stabilizer_tableau()) + state = project_stabilized_state(test_paulis, qvmstab.num_qubits, classical_state=[1, 1]) + state_2 = project_stabilized_state(test_paulis, qvmstab.num_qubits) + + assert np.allclose(np.array(state.todense()), np.array(state_2.todense())) + def test_measurement_noncommuting(): """ @@ -324,18 +383,57 @@ def test_bell_state_measurements(): qvmstab = QVM_Stabilizer() results = qvmstab.run(prog, trials=5000) assert np.isclose(np.mean(results), 0.5, rtol=0.1) + assert all([x[0] == x[1] for x in results]) + + +def test_random_stabilizer_circuit(): + """ + Generate a random stabilizer circuit from {CNOT, H, S, MEASURE} + + Compare the outcome to full state evolution + """ + num_qubits = 2 + # for each qubit pick a set of operations + gate_operations = {1: H, 2: S, 3: CNOT} + np.random.seed(42) + num_gates = 1000 # program has 100 gates in it + + prog = Program() + for _ in range(num_gates): + for jj in range(num_qubits): + instruction_idx = np.random.randint(1, 4) + if instruction_idx == 3: + gate = gate_operations[instruction_idx](jj, (jj + 1) % num_qubits) + else: + gate = gate_operations[instruction_idx](jj) + prog.inst(gate) + + qvmstab = QVM_Stabilizer() + wf = qvmstab.wavefunction(prog) + wf = wf.amplitudes.reshape((-1, 1)) + rho_trial = wf.dot(np.conj(wf.T)) + + qvm = QVMConnection(type_trans='wavefunction') + rho, _ = qvm.wavefunction(prog) + rho = rho.amplitudes.reshape((-1, 1)) + rho_true = rho.dot(np.conj(rho.T)) + assert np.allclose(rho_true, rho_trial) + + rho_from_stab = qvmstab.density(prog) + assert np.allclose(rho_from_stab, rho_true) if __name__ == "__main__": - test_initialization() - test_row_sum_sub_algorithm_g_test() - test_stabilizer_to_matrix_conversion() - test_rowsum_phase_accumulator() - test_rowsum_full() - test_simulation_hadamard() - test_simulation_phase() - test_simulation_cnot() - test_measurement_noncommuting() - test_measurement_commuting() - test_measurement_commuting_result_one() - test_bell_state_measurements() + # test_initialization() + # test_row_sum_sub_algorithm_g_test() + # test_stabilizer_to_matrix_conversion() + # test_rowsum_phase_accumulator() + # test_rowsum_full() + # test_simulation_hadamard() + # test_simulation_phase() + # test_simulation_cnot() + # test_measurement_noncommuting() + # test_measurement_commuting() + # test_measurement_commuting_result_one() + # test_bell_state_measurements() + test_random_stabilizer_circuit() diff --git a/referenceqvm/tests/test_state_actions.py b/referenceqvm/tests/test_state_actions.py index 9240512..ee41623 100644 --- a/referenceqvm/tests/test_state_actions.py +++ b/referenceqvm/tests/test_state_actions.py @@ -260,9 +260,9 @@ def test_stabilizer_projection_ZZZ(): test if we project out the correct state """ stabilizer_state = project_stabilized_state([sZ(0) * sZ(1), sZ(1) * sZ(2), - sX(0) * sX(1) * sX(2)]) + sX(0) * sX(1) * sX(2)]) true_state = np.zeros((8, 1)) true_state[0, 0] = true_state[7, 0] = 1 true_state /= np.sqrt(2) - assert np.allclose(true_state, stabilizer_state.todense()) + assert np.allclose(true_state, np.array(stabilizer_state.todense())) From 93fe02033b3709783fcd35ed7ebd41707c4d13fa Mon Sep 17 00:00:00 2001 From: ncrubin Date: Thu, 2 Aug 2018 10:54:17 -0700 Subject: [PATCH 06/10] Tested CHS simulator returns wavefunction and density matrix. --- referenceqvm/qvm_stabilizer.py | 286 +++--------------- .../{state_actions.py => stabilizer_utils.py} | 87 +++++- ...abilzier_qvm.py => test_stabilizer_qvm.py} | 47 +-- ...te_actions.py => test_stabilizer_utils.py} | 32 +- 4 files changed, 167 insertions(+), 285 deletions(-) rename referenceqvm/{state_actions.py => stabilizer_utils.py} (57%) rename referenceqvm/tests/{test_stabilzier_qvm.py => test_stabilizer_qvm.py} (88%) rename referenceqvm/tests/{test_state_actions.py => test_stabilizer_utils.py} (85%) diff --git a/referenceqvm/qvm_stabilizer.py b/referenceqvm/qvm_stabilizer.py index 00bbd1b..8603091 100644 --- a/referenceqvm/qvm_stabilizer.py +++ b/referenceqvm/qvm_stabilizer.py @@ -15,12 +15,9 @@ # limitations under the License. ############################################################################## """ -Pure QVM that only executes pyQuil programs containing Clifford group generators, -and return the wavefunction or stabilizer. This is based off of the paper by -Aaronson and Gottesman Phys. Rev. A 70, 052328 +Pure QVM that only executes pyQuil programs containing Clifford group generators. +This is based off of the paper by Aaronson and Gottesman Phys. Rev. A 70, 052328 """ -import sys -import numpy as np from functools import reduce from pyquil.quil import Program, get_classical_addresses_from_program from pyquil.quilbase import * @@ -30,7 +27,9 @@ from referenceqvm.qam import QAM from referenceqvm.gates import stabilizer_gate_matrix from referenceqvm.unitary_generator import value_get, tensor_up -from referenceqvm.state_actions import project_stabilized_state +from referenceqvm.stabilizer_utils import (project_stabilized_state, + symplectic_inner_product, + binary_stabilizer_to_pauli_stabilizer) class QVM_Stabilizer(QAM): @@ -79,7 +78,7 @@ def load_program(self, pyquil_program): unitary check QVM checks makes it annoying to allow measurements. TODO: Change load_program for each qvm-subclass - :param Program pyquil_program: program to be ran + :param Program pyquil_program: program to be run """ # typecheck if not isinstance(pyquil_program, Program): @@ -127,8 +126,10 @@ def _n_qubit_tableau(self, num_qubits): """ Construct an empty tableau for a n-qubit system - :param num_qubits: - :return: + :param Int num_qubits: number of qubits represented by the tableau + :return: a numpy integer array representing the tableau stabilizing + the |000..000> state + :rtype: np.ndarray """ tableau = np.zeros((2 * num_qubits, (2 * num_qubits) + 1), dtype=int) @@ -143,7 +144,6 @@ def transition(self, instruction): Implements a full transition, including pre/post noise hooks. :param QuilAction instruction: instruction to be executed - :return: if program is halted TRUE, else FALSE :rtype: bool int """ @@ -203,15 +203,17 @@ def _transition(self, instruction): def run(self, pyquil_program, classical_addresses=None, trials=1): """ - Run program. + Run pyquil program and return the results Loads and checks program if all gates are within the stabilizer set. - Then executes program + Then executes program. If measurements are requested then the measured + results are returned - :param porgram: - :param classical_addresses: - :param trials: - :return: + :param Program pyquil_program: a pyquil Program containing only + CNOT-H-S-MEASUREMENT operations + :param classical_addresses: classical addresses to return + :param trials: number of times to repeat the execution of the program + :return: list of lists of classical memory after each run """ self.load_program(pyquil_program) @@ -236,9 +238,11 @@ def density(self, pyquil_program): Run program and compute the density matrix Loads and checks program if all gates are within the stabilizer set. - Then executes program + Then executes program and returns the final density matrix for the + stabilizer - :param porgram: + :param Program pyquil_program: a pyquil Program containing only + CNOT-H-S-MEASUREMENT operations :return: """ self.load_program(pyquil_program) @@ -251,7 +255,7 @@ def density(self, pyquil_program): def wavefunction(self, program): """ - Simulate the stabilizer program and then project out the final state. + Simulate the program and then project out the final state. Return the final state as a wavefunction pyquil object @@ -288,7 +292,7 @@ def _rowsum(self, h, i): :param Int h: row index 'h' :param Int i: row index 'i' - :return: + :return: None but mutates the tableau """ # NOTE: this is left multiplication of P(i). P(i) * P(h) phase_accumulator = self._rowsum_phase_accumulator(h, i) @@ -311,6 +315,7 @@ def _rowsum_phase_accumulator(self, h, i): phase accumulator sub algorithm for the rowsum routine Note: this accumulates the $i$ phase for row_i * row_h NOT row_h, row_i + :param Int h: row index 'h' :param Int i: row index 'i' :return: phase mod 4 @@ -332,11 +337,11 @@ def _g_update(self, x1, z1, x2, z2): when the pauli matrices represented by x1z1 and x2z2 are multiplied together - :param x1: - :param z1: - :param x2: - :param z2: - :return: + :param x1: binary variable for the first pauli operator representing X + :param z1: binary variable for th first pauli operator representing Z + :param x2: binary variable for the second pauli operator representing X + :param z2: binary variable for the second pauli operator representing Z + :return: power that 1j is raised to when multiplying paulis together """ # if the first term is identity if x1 == z1 == 0: @@ -366,13 +371,14 @@ def _g_update(self, x1, z1, x2, z2): if x1 == 0 and z1 == 1: return x2 * (1 - 2 * z2) - raise ValueError("we were unable to multiply the pauli operators together!") + raise ValueError("We were unable to multiply the pauli operators" + + " together!") def _apply_cnot(self, instruction): """ Apply a CNOT to the tableau - :param instruction: pyquil abstract instruction. Must have + :param instruction: pyquil abstract instruction. """ a, b = [value_get(x) for x in instruction.qubits] for i in range(2 * self.num_qubits): @@ -382,7 +388,7 @@ def _apply_cnot(self, instruction): def _cnot_phase_update(self, i, c, t): """ - update r_{i} + update r_{i} as a submethod to applying CNOT to the tableau :param i: tableau row index :param c: control qubit index @@ -396,8 +402,7 @@ def _apply_hadamard(self, instruction): """ Apply a hadamard gate on qubit defined in instruction - :param instruction: - :return: + :param instruction: pyquil abstract instruction. """ qubit_label = [value_get(x) for x in instruction.qubits][0] for i in range(2 * self.num_qubits): @@ -408,8 +413,7 @@ def _apply_phase(self, instruction): """ Apply the phase gate instruction ot the tableau - :param instruction: - :return: + :param instruction: pyquil abstract instruction. """ qubit_label = [value_get(x) for x in instruction.qubits][0] for i in range(2 * self.num_qubits): @@ -417,6 +421,11 @@ def _apply_phase(self, instruction): self.tableau[i, qubit_label + self.num_qubits] = self.tableau[i, qubit_label + self.num_qubits] ^ self.tableau[i, qubit_label] def _apply_measurement(self, instruction): + """ + Apply a measurement + + :param instruction: pyquil abstract instruction. + """ t_qbit = value_get(instruction.qubit) t_cbit = value_get(instruction.classical_reg) @@ -424,9 +433,9 @@ def _apply_measurement(self, instruction): # this is analogous to the case when the measurement operator does not # commute with at least one stabilizer if any(self.tableau[self.num_qubits:, t_qbit] == 1): - # find the first one. + # find the first `1'. xpa_idx = np.where(self.tableau[self.num_qubits:, t_qbit] == 1)[0][0] # take the first index - xpa_idx += self.num_qubits # adjust so we can index into the tableau + xpa_idx += self.num_qubits # adjust so we can index into the tableau for ii in range(2 * self.num_qubits): # loop over each row and call rowsum(ii, xpa_idx) if ii != xpa_idx and self.tableau[ii, t_qbit] == 1: self._rowsum(ii, xpa_idx) @@ -447,223 +456,24 @@ def _apply_measurement(self, instruction): # outcome of measurement is deterministic...need to determine sign else: - # augment tableaue with a scratch space + # augment tableau with a scratch space self.tableau = np.vstack((self.tableau, np.zeros((1, 2 * self.num_qubits + 1), dtype=int))) for ii in range(self.num_qubits): # We need to check if R(i) anticommutes with Za...which it does if x_{ia} = 1 - if self.tableau[ii, t_qbit] == 1: # refrencing the destabilizers + if self.tableau[ii, t_qbit] == 1: # referencing the destabilizers + # TODO: Remove this for performance? # check something silly. Does the destabilizer anticommute with the observable? It SHOULD! tmp_vector_representing_z_qubit = np.zeros((2 * self.num_qubits), dtype=int) tmp_vector_representing_z_qubit[t_qbit] = 1 assert symplectic_inner_product(tmp_vector_representing_z_qubit, self.tableau[ii, :-1]) == 1 # row sum on the stabilizers (summing up operators such that we get operator Z_{a}) - self._rowsum(2 * self.num_qubits, ii + self.num_qubits) # note: A-G says 2 n + 1...this is correct...but they start counting at 1 not zero + # note: A-G says 2 n + 1...this is correct...but they start counting at 1 not zero + self._rowsum(2 * self.num_qubits, ii + self.num_qubits) # set the classical bit to be the last element of the scratch row self.classical_memory[t_cbit] = self.tableau[-1, -1] # remove the scratch space self.tableau = self.tableau[:2 * self.num_qubits, :] - - - -def pauli_stabilizer_to_binary_stabilizer(stabilizer_list): - """ - Convert a list of stabilizers represented as PauliTerms to a binary tableau form - - :param List stabilizer_list: list of stabilizers where each element is a PauliTerm - :return: return an integer matrix representing the stabilizers where each row is a - stabilizer. The size of the matrix is n x (2 * n) where n is the maximum - qubit index. - """ - if not all([isinstance(x, PauliTerm) for x in stabilizer_list]): - raise TypeError("At least one element of stabilizer_list is not a PauliTerm") - - max_qubit = max([max(x.get_qubits()) for x in stabilizer_list]) + 1 - stabilizer_tableau = np.zeros((len(stabilizer_list), 2 * max_qubit + 1), dtype=int) - for row_idx, term in enumerate(stabilizer_list): - for i, pterm in term: # iterator for each tensor-product element of the Pauli operator - if pterm == 'X': - stabilizer_tableau[row_idx, i] = 1 - elif pterm == 'Z': - stabilizer_tableau[row_idx, i + max_qubit] = 1 - elif pterm == 'Y': - stabilizer_tableau[row_idx, i] = 1 - stabilizer_tableau[row_idx, i + max_qubit] = 1 - else: - # term is identity - pass - - if not (np.isclose(term.coefficient, -1) or np.isclose(term.coefficient, 1)): - raise ValueError("stabilizers must have a +/- coefficient") - - if int(np.sign(term.coefficient.real)) == 1: - stabilizer_tableau[row_idx, -1] = 0 - elif int(np.sign(term.coefficient.real)) == -1: - stabilizer_tableau[row_idx, -1] = 1 - else: - raise TypeError('unrecognized on pauli term of stabilizer') - - return stabilizer_tableau - - -def binary_stabilizer_to_pauli_stabilizer(stabilizer_tableau): - """ - Convert a stabilizer tableau to a list of PauliTerms - - :param stabilizer_tableau: - :return: - """ - stabilizer_list = [] - num_qubits = (stabilizer_tableau.shape[1] - 1) // 2 - for nn in range(stabilizer_tableau.shape[0]): # iterate through the rows - stabilizer_element = [] - for ii in range(num_qubits): - if stabilizer_tableau[nn, ii] == 1 and stabilizer_tableau[nn, ii + num_qubits] == 0: - stabilizer_element.append(sX(ii)) - elif stabilizer_tableau[nn, ii] == 0 and stabilizer_tableau[nn, ii + num_qubits] == 1: - stabilizer_element.append(sZ(ii)) - elif stabilizer_tableau[nn, ii] == 1 and stabilizer_tableau[nn, ii + num_qubits] == 1: - stabilizer_element.append(sY(ii)) - - stabilizer_term = reduce(lambda x, y: x * y, stabilizer_element) * ((-1) ** stabilizer_tableau[nn, -1]) - stabilizer_list.append(stabilizer_term) - return stabilizer_list - - -def binary_rref(code_matrix): - """ - Convert the binary code_matrix into rref - - :param code_matrix: - :return: - """ - code_matrix = code_matrix.astype(int) - rdim, cdim = code_matrix.shape - for ridx in range(rdim): - # print("start [{},{}] = {}".format(ridx, ridx, - # code_matrix[ridx, ridx])) - # set first element - if not np.isclose(code_matrix[ridx, ridx], 1): - for ii in range(ridx + 1, rdim): - if np.isclose(code_matrix[ii, ridx], 1.0): - # switch rows - # code_matrix[ridx, :], code_matrix[ii, :] = code_matrix[ii, :], code_matrix[ridx, :] - code_matrix[[ridx, ii]] = code_matrix[[ii, ridx]] - print("r_{} <-> r_{}".format(ii + 1, ridx + 1)) - - # print(code_matrix) - # start elimination of all other rows - # print("pivot [{}, {}] = {}".format(ridx, ridx, - # code_matrix[ridx, ridx])) - - rows_to_eliminate = set(range(rdim)) - rows_to_eliminate.remove(ridx) - # print(rows_to_eliminate) - for elim_ridx in list(rows_to_eliminate): - # eliminate by Hadamard product on the rows (XOR element wise) - # if other row element is 1 then eliminate by subtracting ridx row - - # if np.isclose(code_matrix[elim_ridx, ridx], 1.0): - # code_matrix[elim_ridx, :] = code_matrix[elim_ridx, :] - \ - # code_matrix[ridx, :] - if np.isclose(code_matrix[elim_ridx, ridx], 1.0): - code_matrix[elim_ridx, :] = code_matrix[elim_ridx, :] ^ code_matrix[ridx, :] - # print("elim [{}, :]".format(elim_ridx)) - - # print("semi rref at row {}".format(ridx)) - # print(code_matrix) - return code_matrix - - -def symplectic_inner_product(vector1, vector2): - """ - Operators commute if their binary form symplectic inner product is zero - - Operators anticommute if their binary form symplectic inner product is one - - :param vector1: binary form of operator with no sign info - :param vector2: binary form of a pauli operator with no sign info - :return: 0, 1 - """ - if vector1.shape != vector2.shape: - raise ValueError("vectors must be the same size.") - - # TODO: add a check for binary or integer linear arrays - - hadamard_product = np.multiply(vector1, vector2) - return reduce(lambda x, y: x ^ y, hadamard_product) - - -if __name__ == "__main__": - # practice my clifford transformations - # = |0000...0N> - # N |psi> = N g |psi> = NgNd N|psi> = g2 |psi2> = |psi2> - - from referenceqvm.gates import X, Y, Z, H, CNOT, S, I - from sympy import Matrix - - # HZH = X - test_X = H.dot(Z).dot(H) - assert np.allclose(test_X, X) - - # HXH = Z - test_Z = H.dot(X).dot(H) - assert np.allclose(test_Z, Z) - - # HYH = iHZHHXH = iXZ = -Y - test_nY = H.dot(Y).dot(H) - assert np.allclose(test_nY, -Y) - - # XYX = -Y - assert np.allclose(X.dot(Y).dot(X), -Y) - - # XZX = -Z - assert np.allclose(X.dot(Z).dot(X), -Z) - - # SXS = Y - assert np.allclose(S.dot(X).dot(np.conj(S).T), Y) - - # CNOT(X otimes X)CNOT = X otimes I - # (C, T) for CNOT as written in reference-qvm - assert np.allclose(CNOT.dot(np.kron(I, X)).dot(CNOT), np.kron(I, X)) - assert np.allclose(CNOT.dot(np.kron(X, I)).dot(CNOT), np.kron(X, X)) - assert np.allclose(CNOT.dot(np.kron(X, X)).dot(CNOT), np.kron(X, I)) - - # CNOT(Z otimes Z)CNOT = Z otimes I - assert np.allclose(CNOT.dot(np.kron(Z, Z)).dot(CNOT), np.kron(I, Z)) - assert np.allclose(CNOT.dot(np.kron(Z, I)).dot(CNOT), np.kron(Z, I)) - assert np.allclose(CNOT.dot(np.kron(I, Z)).dot(CNOT), np.kron(Z, Z)) - - Z1 = np.kron(Z, I) - Z2 = np.kron(I, Z) - ZZ = np.kron(Z, Z) - X1 = np.kron(X, I) - X2 = np.kron(I, X) - XX = np.kron(X, X) - YY = np.kron(Y, Y) - - generator_matrix = np.zeros((3, 4)) - generator_matrix[0, 2] = generator_matrix[0, 3] = 1 # ZZ - generator_matrix[1, 0] = generator_matrix[1, 1] = 1 # XX - generator_matrix[2, 0] = generator_matrix[2, 2] = 1 # Y - generator_matrix[2, 1] = generator_matrix[2, 3] = 1 # Y - parity_vector = np.zeros((3, 1)) - parity_vector[-1] = 1 - - # define the amplitudes from the generator set by diagonalizing - # the tableau matrix? - generator_matrix = np.hstack((generator_matrix, parity_vector)) - print(generator_matrix) - - generator_matrix_2 = binary_rref(generator_matrix) - print('\n') - print(generator_matrix_2) - - mat = Matrix(generator_matrix) - mat, pivots = mat.rref() - mat = np.asarray(mat).astype(int) - print(mat) - diff --git a/referenceqvm/state_actions.py b/referenceqvm/stabilizer_utils.py similarity index 57% rename from referenceqvm/state_actions.py rename to referenceqvm/stabilizer_utils.py index 4253875..136c290 100644 --- a/referenceqvm/state_actions.py +++ b/referenceqvm/stabilizer_utils.py @@ -10,9 +10,10 @@ Given """ import sys +from functools import reduce from scipy.sparse import csc_matrix import numpy as np -from pyquil.paulis import PauliTerm +from pyquil.paulis import PauliTerm, sX, sZ, sY def compute_action(classical_state, pauli_operator, num_qubits): @@ -140,3 +141,87 @@ def project_stabilized_state(stabilizer_list, num_qubits=None, state /= np.sqrt(float(normalization.real)) # this is needed or it will cast as a np.matrix return state + + +def pauli_stabilizer_to_binary_stabilizer(stabilizer_list): + """ + Convert a list of stabilizers represented as PauliTerms to a binary tableau form + + :param List stabilizer_list: list of stabilizers where each element is a PauliTerm + :return: return an integer matrix representing the stabilizers where each row is a + stabilizer. The size of the matrix is n x (2 * n) where n is the maximum + qubit index. + """ + if not all([isinstance(x, PauliTerm) for x in stabilizer_list]): + raise TypeError("At least one element of stabilizer_list is not a PauliTerm") + + max_qubit = max([max(x.get_qubits()) for x in stabilizer_list]) + 1 + stabilizer_tableau = np.zeros((len(stabilizer_list), 2 * max_qubit + 1), dtype=int) + for row_idx, term in enumerate(stabilizer_list): + for i, pterm in term: # iterator for each tensor-product element of the Pauli operator + if pterm == 'X': + stabilizer_tableau[row_idx, i] = 1 + elif pterm == 'Z': + stabilizer_tableau[row_idx, i + max_qubit] = 1 + elif pterm == 'Y': + stabilizer_tableau[row_idx, i] = 1 + stabilizer_tableau[row_idx, i + max_qubit] = 1 + else: + # term is identity + pass + + if not (np.isclose(term.coefficient, -1) or np.isclose(term.coefficient, 1)): + raise ValueError("stabilizers must have a +/- coefficient") + + if int(np.sign(term.coefficient.real)) == 1: + stabilizer_tableau[row_idx, -1] = 0 + elif int(np.sign(term.coefficient.real)) == -1: + stabilizer_tableau[row_idx, -1] = 1 + else: + raise TypeError('unrecognized on pauli term of stabilizer') + + return stabilizer_tableau + + +def binary_stabilizer_to_pauli_stabilizer(stabilizer_tableau): + """ + Convert a stabilizer tableau to a list of PauliTerms + + :param stabilizer_tableau: Stabilizer tableau to turn into pauli terms + :return: a list of PauliTerms representing the tableau + :rytpe: List of PauliTerms + """ + stabilizer_list = [] + num_qubits = (stabilizer_tableau.shape[1] - 1) // 2 + for nn in range(stabilizer_tableau.shape[0]): # iterate through the rows + stabilizer_element = [] + for ii in range(num_qubits): + if stabilizer_tableau[nn, ii] == 1 and stabilizer_tableau[nn, ii + num_qubits] == 0: + stabilizer_element.append(sX(ii)) + elif stabilizer_tableau[nn, ii] == 0 and stabilizer_tableau[nn, ii + num_qubits] == 1: + stabilizer_element.append(sZ(ii)) + elif stabilizer_tableau[nn, ii] == 1 and stabilizer_tableau[nn, ii + num_qubits] == 1: + stabilizer_element.append(sY(ii)) + + stabilizer_term = reduce(lambda x, y: x * y, stabilizer_element) * ((-1) ** stabilizer_tableau[nn, -1]) + stabilizer_list.append(stabilizer_term) + return stabilizer_list + + +def symplectic_inner_product(vector1, vector2): + """ + Operators commute if the symplectic inner product of their binary form is zero + + Operators anticommute if symplectic inner product of their binary form is one + + :param vector1: binary form of operator with no sign info + :param vector2: binary form of a pauli operator with no sign info + :return: 0, 1 + """ + if vector1.shape != vector2.shape: + raise ValueError("vectors must be the same size.") + + # TODO: add a check for binary or integer linear arrays + + hadamard_product = np.multiply(vector1, vector2) + return reduce(lambda x, y: x ^ y, hadamard_product) diff --git a/referenceqvm/tests/test_stabilzier_qvm.py b/referenceqvm/tests/test_stabilizer_qvm.py similarity index 88% rename from referenceqvm/tests/test_stabilzier_qvm.py rename to referenceqvm/tests/test_stabilizer_qvm.py index c915315..fbece46 100644 --- a/referenceqvm/tests/test_stabilzier_qvm.py +++ b/referenceqvm/tests/test_stabilizer_qvm.py @@ -8,9 +8,10 @@ from pyquil.paulis import sI, sX, sY, sZ, PAULI_COEFF, PAULI_OPS from pyquil.quil import Program from pyquil.gates import H, S, CNOT, I -from referenceqvm.qvm_stabilizer import (QVM_Stabilizer, pauli_stabilizer_to_binary_stabilizer, - binary_stabilizer_to_pauli_stabilizer) -from referenceqvm.state_actions import project_stabilized_state +from referenceqvm.qvm_stabilizer import QVM_Stabilizer +from referenceqvm.stabilizer_utils import (pauli_stabilizer_to_binary_stabilizer, + binary_stabilizer_to_pauli_stabilizer) +from referenceqvm.stabilizer_utils import project_stabilized_state from referenceqvm.api import QVMConnection @@ -43,30 +44,6 @@ def test_initialization(): assert np.allclose(initial_tableau, qvmstab.tableau) -def test_stabilizer_to_matrix_conversion(): - # bitflip code - stabilizer_matrix = pauli_stabilizer_to_binary_stabilizer(bell_stabilizer) - true_stabilizer_matrix = np.array([[0, 0, 1, 1, 0], - [1, 1, 0, 0, 0]]) - assert np.allclose(true_stabilizer_matrix, stabilizer_matrix) - - test_stabilizer_list = binary_stabilizer_to_pauli_stabilizer(true_stabilizer_matrix) - for idx, term in enumerate(test_stabilizer_list): - assert term == bell_stabilizer[idx] - - # given some codes convert them to code matrices - stabilizer_matrix = pauli_stabilizer_to_binary_stabilizer(five_qubit_code_generators) - true_stabilizer_matrix = np.array([[1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0], - [0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0], - [1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0], - [0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0]]) - assert np.allclose(true_stabilizer_matrix, stabilizer_matrix) - - test_stabilizer_list = binary_stabilizer_to_pauli_stabilizer(true_stabilizer_matrix) - for idx, term in enumerate(test_stabilizer_list): - assert term == five_qubit_code_generators[idx] - - def test_row_sum_sub_algorithm_g_test(): """ Test the behavior of the subroutine of the rowsum routine @@ -421,19 +398,3 @@ def test_random_stabilizer_circuit(): rho_from_stab = qvmstab.density(prog) assert np.allclose(rho_from_stab, rho_true) - - -if __name__ == "__main__": - # test_initialization() - # test_row_sum_sub_algorithm_g_test() - # test_stabilizer_to_matrix_conversion() - # test_rowsum_phase_accumulator() - # test_rowsum_full() - # test_simulation_hadamard() - # test_simulation_phase() - # test_simulation_cnot() - # test_measurement_noncommuting() - # test_measurement_commuting() - # test_measurement_commuting_result_one() - # test_bell_state_measurements() - test_random_stabilizer_circuit() diff --git a/referenceqvm/tests/test_state_actions.py b/referenceqvm/tests/test_stabilizer_utils.py similarity index 85% rename from referenceqvm/tests/test_state_actions.py rename to referenceqvm/tests/test_stabilizer_utils.py index ee41623..039aebe 100644 --- a/referenceqvm/tests/test_state_actions.py +++ b/referenceqvm/tests/test_stabilizer_utils.py @@ -3,8 +3,11 @@ eigenspace of a set of generators or stabilizers """ import numpy as np -from referenceqvm.state_actions import compute_action, project_stabilized_state -from referenceqvm.tests.test_stabilzier_qvm import five_qubit_code_generators +from referenceqvm.stabilizer_utils import (compute_action, project_stabilized_state, + binary_stabilizer_to_pauli_stabilizer, + pauli_stabilizer_to_binary_stabilizer) +from referenceqvm.tests.test_stabilizer_qvm import (five_qubit_code_generators, + bell_stabilizer) from pyquil.paulis import sX, sZ, sY, sI, PauliSum import pytest @@ -29,6 +32,30 @@ def test_compute_action_type_checks(): compute_action('0001', sX(0), 4) +def test_stabilizer_to_matrix_conversion(): + # bitflip code + stabilizer_matrix = pauli_stabilizer_to_binary_stabilizer(bell_stabilizer) + true_stabilizer_matrix = np.array([[0, 0, 1, 1, 0], + [1, 1, 0, 0, 0]]) + assert np.allclose(true_stabilizer_matrix, stabilizer_matrix) + + test_stabilizer_list = binary_stabilizer_to_pauli_stabilizer(true_stabilizer_matrix) + for idx, term in enumerate(test_stabilizer_list): + assert term == bell_stabilizer[idx] + + # given some codes convert them to code matrices + stabilizer_matrix = pauli_stabilizer_to_binary_stabilizer(five_qubit_code_generators) + true_stabilizer_matrix = np.array([[1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0], + [0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0], + [1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0], + [0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0]]) + assert np.allclose(true_stabilizer_matrix, stabilizer_matrix) + + test_stabilizer_list = binary_stabilizer_to_pauli_stabilizer(true_stabilizer_matrix) + for idx, term in enumerate(test_stabilizer_list): + assert term == five_qubit_code_generators[idx] + + def test_compute_action_identity(): """ Action of Pauli operators on state @@ -265,4 +292,3 @@ def test_stabilizer_projection_ZZZ(): true_state[0, 0] = true_state[7, 0] = 1 true_state /= np.sqrt(2) assert np.allclose(true_state, np.array(stabilizer_state.todense())) - From 34c31eb61946abb3a5c31cf34b77e70888d790d2 Mon Sep 17 00:00:00 2001 From: ncrubin Date: Thu, 2 Aug 2018 16:32:49 -0700 Subject: [PATCH 07/10] update doc strings --- referenceqvm/stabilizer_utils.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/referenceqvm/stabilizer_utils.py b/referenceqvm/stabilizer_utils.py index 136c290..58333de 100644 --- a/referenceqvm/stabilizer_utils.py +++ b/referenceqvm/stabilizer_utils.py @@ -9,7 +9,6 @@ Given """ -import sys from functools import reduce from scipy.sparse import csc_matrix import numpy as np @@ -73,9 +72,10 @@ def state_family_generator(state, pauli_operator): This is accomplished in a sparse format where a sparse vector is returned after the action is accumulate in a new list of data and indices - :param state: - :param pauli_operator: - :return: + :param csc_matrix state: wavefunction represented as a column sparse matrix + :param PauliTerm pauli_operator: action to apply to the state + :return: new state + :rtype: csc_matrix """ if not isinstance(state, csc_matrix): raise TypeError("we only take csc_matrix") @@ -104,10 +104,9 @@ def project_stabilized_state(stabilizer_list, num_qubits=None, |psi> = (1/2^{n}) * Product_{i=0}{n-1}[ 1 + G_{i}] |vac> - :param stabilizer_list: + :param List stabilizer_list: set of PauliTerms that are the stabilizers :param num_qubits: integer number of qubits - :param classical_state: Default None. Defaults to |+>^{otimes n} - + :param classical_state: Default None. Defaults to |+>^{\otimes n} :return: state projected by stabilizers """ if num_qubits is None: From c66b30f28272c0dc8e338864acdd5aa354d39b3c Mon Sep 17 00:00:00 2001 From: ncrubin Date: Thu, 2 Aug 2018 16:35:45 -0700 Subject: [PATCH 08/10] Fixed imports --- referenceqvm/api.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/referenceqvm/api.py b/referenceqvm/api.py index 5d6d4db..b7e2991 100644 --- a/referenceqvm/api.py +++ b/referenceqvm/api.py @@ -1,10 +1,11 @@ """ Sets up the appropriate QVM, and interfacing. """ -from referenceqvm.gates import gate_matrix +from referenceqvm.gates import gate_matrix, stabilizer_gate_matrix from referenceqvm.qvm_wavefunction import QVM_Wavefunction from referenceqvm.qvm_unitary import QVM_Unitary from referenceqvm.qvm_density import QVM_Density +from referenceqvm.qvm_stabilizer import QVM_Stabilizer def QVMConnection(type_trans='wavefunction', @@ -36,7 +37,7 @@ def QVMConnection(type_trans='wavefunction', qvm = QVM_Density(gate_set=gate_set, noise_model=noise_model) elif type_trans == 'stabilizer': - qvm = QVM_Stabilizer(gate_set=stabilizer_gate_set) + qvm = QVM_Stabilizer(gate_set=stabilizer_gate_matrix) else: raise TypeError("{} is not a valid QVM type.".format(type_trans)) From 293665db13a9aa23bdf326f96a9b867222d34691 Mon Sep 17 00:00:00 2001 From: ncrubin Date: Thu, 2 Aug 2018 16:37:52 -0700 Subject: [PATCH 09/10] Remove slot object check Slots are no longer in pyquil --- referenceqvm/unitary_generator.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/referenceqvm/unitary_generator.py b/referenceqvm/unitary_generator.py index 90e1733..ce2f9ce 100644 --- a/referenceqvm/unitary_generator.py +++ b/referenceqvm/unitary_generator.py @@ -21,7 +21,6 @@ Note: uses SciPy sparse diagonal (DIA) representation to increase space and timeefficiency. """ -import numpy as np from collections import Sequence from numbers import Integral @@ -418,8 +417,6 @@ def value_get(param_obj): return param_obj.index elif isinstance(param_obj, Addr): return param_obj.address - # elif isinstance(param_obj, Slot): - # return param_obj.value() elif isinstance(param_obj, Label): return param_obj.name else: From 6a25f86787cbab40ab6666d91d1899b3bfba1b64 Mon Sep 17 00:00:00 2001 From: ncrubin Date: Thu, 9 Aug 2018 12:13:57 -0700 Subject: [PATCH 10/10] PR updates and improvements Clean up and list comps --- referenceqvm/qam.py | 23 ++++++++++++++++++++++- referenceqvm/qvm_stabilizer.py | 26 +++++++++++++------------- referenceqvm/qvm_wavefunction.py | 20 -------------------- referenceqvm/stabilizer_utils.py | 6 +++--- 4 files changed, 38 insertions(+), 37 deletions(-) diff --git a/referenceqvm/qam.py b/referenceqvm/qam.py index 05a71ee..27cad76 100644 --- a/referenceqvm/qam.py +++ b/referenceqvm/qam.py @@ -27,7 +27,8 @@ from pyquil.quilbase import (Gate, Measurement, UnaryClassicalInstruction, - BinaryClassicalInstruction) + BinaryClassicalInstruction, + Label, JumpTarget) class QAM(object): @@ -179,6 +180,26 @@ def kernel(self): if halted: break + def find_label(self, label): + """ + Helper function that iterates over the program and looks for a + JumpTarget that has a Label matching the input label. + + :param Label label: Label object to search for in program + + :return: program index where Label is found + :rtype: int + """ + assert isinstance(label, Label) + for index, action in enumerate(self.program): + if isinstance(action, JumpTarget): + if label == action.label: + return index + + # Label was not found in program. + raise RuntimeError("Improper program - Jump Target not found in the " + "input program!") + def transition(self, instruction): """ Abstract class for the transition type diff --git a/referenceqvm/qvm_stabilizer.py b/referenceqvm/qvm_stabilizer.py index 8603091..d619e13 100644 --- a/referenceqvm/qvm_stabilizer.py +++ b/referenceqvm/qvm_stabilizer.py @@ -88,17 +88,16 @@ def load_program(self, pyquil_program): raise NotImplementedError("QAM needs to be subclassed in order to " "load program") - # create defgate dictionary - defined_gates = {} - for dg in pyquil_program.defined_gates: - defined_gates[dg.name] = dg.matrix - self.defgate_set = defined_gates + # # create defgate dictionary + # defined_gates = {} + # for dg in pyquil_program.defined_gates: + # defined_gates[dg.name] = dg.matrix + # self.defgate_set = defined_gates - # if QVM_Unitary, check if all instructions are valid. invalid = False for instr in pyquil_program: if isinstance(instr, Gate): - if not (instr.name in self.gate_set.keys() or instr.name in self.defgate_set.keys()): + if not instr.name in self.gate_set.keys(): invalid = True break elif isinstance(instr, Measurement): @@ -226,6 +225,8 @@ def run(self, pyquil_program, classical_addresses=None, trials=1): self.tableau = self._n_qubit_tableau(self.num_qubits) self.kernel() results.append(list(map(int, self.classical_memory[classical_addresses]))) + # results.append([int(b) for b in self.classical_memory[classical_addresses]]) + assert results[-1] == [int(b) for b in self.classical_memory[classical_addresses]] # reset qvm self.memory_reset() @@ -249,8 +250,7 @@ def density(self, pyquil_program): self.tableau = self._n_qubit_tableau(self.num_qubits) self.kernel() stabilizers = binary_stabilizer_to_pauli_stabilizer(self.stabilizer_tableau()) - pauli_ops = list(map(lambda x: 0.5 * (sI(0) + x), stabilizers)) - pauli_ops = reduce(lambda x, y: x * y, pauli_ops) + pauli_ops = reduce(lambda x, y: x * y, [0.5 * (sI(0) + b) for b in stabilizers]) return tensor_up(pauli_ops, self.num_qubits) def wavefunction(self, program): @@ -267,7 +267,7 @@ def wavefunction(self, program): self.kernel() stabilizers = binary_stabilizer_to_pauli_stabilizer(self.stabilizer_tableau()) stabilizer_state = project_stabilized_state(stabilizers) - stabilizer_state = np.array(stabilizer_state.todense()) # todense() returns a matrixlib type + stabilizer_state = stabilizer_state.toarray() return Wavefunction(stabilizer_state.flatten()) def stabilizer_tableau(self): @@ -322,7 +322,7 @@ def _rowsum_phase_accumulator(self, h, i): """ phase_accumulator = 0 for j in range(self.num_qubits): - # self._g_update(x_{hj}, z_{hj}, x_{ij}, z_{ij}) + # self._g_update(x_{ij}, z_{ij}, x_{hj}, z_{hj}) phase_accumulator += self._g_update(self.tableau[i, j], self.tableau[i, self.num_qubits + j], self.tableau[h, j], @@ -395,8 +395,8 @@ def _cnot_phase_update(self, i, c, t): :param t: target qubit index :return: 0/1 phase update for r_{i} """ - return self.tableau[i, -1] ^ (self.tableau[i, c] * self.tableau[i, t + self.num_qubits]) * ( - self.tableau[i, t] ^ self.tableau[i, c + self.num_qubits] ^ 1) + return self.tableau[i, -1] ^ (self.tableau[i, c] * self.tableau[i, t + self.num_qubits]) * \ + (self.tableau[i, t] ^ self.tableau[i, c + self.num_qubits] ^ 1) def _apply_hadamard(self, instruction): """ diff --git a/referenceqvm/qvm_wavefunction.py b/referenceqvm/qvm_wavefunction.py index b205927..cc733a1 100644 --- a/referenceqvm/qvm_wavefunction.py +++ b/referenceqvm/qvm_wavefunction.py @@ -128,26 +128,6 @@ def measurement(self, qubit_index, psi=None): return measured_val, unitary - def find_label(self, label): - """ - Helper function that iterates over the program and looks for a - JumpTarget that has a Label matching the input label. - - :param Label label: Label object to search for in program - - :return: program index where Label is found - :rtype: int - """ - assert isinstance(label, Label) - for index, action in enumerate(self.program): - if isinstance(action, JumpTarget): - if label == action.label: - return index - - # Label was not found in program. - raise RuntimeError("Improper program - Jump Target not found in the " - "input program!") - def _transition(self, instruction): """ Implements a transition on the wf-qvm. diff --git a/referenceqvm/stabilizer_utils.py b/referenceqvm/stabilizer_utils.py index 58333de..750a365 100644 --- a/referenceqvm/stabilizer_utils.py +++ b/referenceqvm/stabilizer_utils.py @@ -88,7 +88,7 @@ def state_family_generator(state, pauli_operator): rindices, cindices = state.nonzero() for ridx, cidx in zip(rindices, cindices): # this is so gross looking - bitstring = list(map(int, np.binary_repr(ridx, width=num_qubits)))[::-1] + bitstring = [int(x) for x in np.binary_repr(ridx, width=num_qubits)][::-1] new_ket, new_coefficient = compute_action(bitstring, pauli_operator, num_qubits) new_indices.append(int("".join([str(x) for x in new_ket[::-1]]), 2)) new_coeffs.append(state[ridx, cidx] * new_coefficient * pauli_operator.coefficient) @@ -136,8 +136,8 @@ def project_stabilized_state(stabilizer_list, num_qubits=None, state += state_family_generator(state, generator) state /= 2 - normalization = (state.conj().T.dot(state)).todense() - state /= np.sqrt(float(normalization.real)) # this is needed or it will cast as a np.matrix + normalization = (state.conj().T.dot(state)).toarray() + state /= np.sqrt(float(normalization.real)) return state