diff --git a/referenceqvm/api.py b/referenceqvm/api.py index 422770c..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', @@ -35,6 +36,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_matrix) 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/qam.py b/referenceqvm/qam.py index c9e3507..27cad76 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 @@ -26,7 +27,8 @@ from pyquil.quilbase import (Gate, Measurement, UnaryClassicalInstruction, - BinaryClassicalInstruction) + BinaryClassicalInstruction, + Label, JumpTarget) class QAM(object): @@ -91,8 +93,7 @@ def load_program(self, pyquil_program): # 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 +106,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 @@ -168,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 new file mode 100644 index 0000000..d619e13 --- /dev/null +++ b/referenceqvm/qvm_stabilizer.py @@ -0,0 +1,479 @@ +#!/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 Clifford group generators. +This is based off of the paper by Aaronson and Gottesman Phys. Rev. A 70, 052328 +""" +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.stabilizer_utils import (project_stabilized_state, + symplectic_inner_product, + binary_stabilizer_to_pauli_stabilizer) + + +class QVM_Stabilizer(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 + + S I M U L A T I N G + + 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=stabilizer_gate_matrix, + defgate_set=None): + """ + Subclassed from QAM this is a pure QVM. + """ + 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) + # 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 run + """ + # 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 + + invalid = False + for instr in pyquil_program: + if isinstance(instr, Gate): + if not instr.name in self.gate_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 + + :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) + + # set up initial state |0>^{otimes n} + for ii in range(2 * self.num_qubits): + tableau[ii, ii] = 1 + return tableau + + 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): + """ + 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 pyquil program and return the results + + Loads and checks program if all gates are within the stabilizer set. + Then executes program. If measurements are requested then the measured + results are returned + + :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) + + 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]))) + # 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() + self.program_counter = 0 + + 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 and returns the final density matrix for the + stabilizer + + :param Program pyquil_program: a pyquil Program containing only + CNOT-H-S-MEASUREMENT operations + :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 * y, [0.5 * (sI(0) + b) for b in stabilizers]) + return tensor_up(pauli_ops, self.num_qubits) + + def wavefunction(self, program): + """ + Simulate the 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 = stabilizer_state.toarray() + return Wavefunction(stabilizer_state.flatten()) + + 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: 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) + # 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 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 _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 + """ + phase_accumulator = 0 + for j in range(self.num_qubits): + # 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], + 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 + + def _g_update(self, x1, z1, x2, z2): + """ + function that takes 4 bits and returns the power of $i$ {0, 1, -1} + + when the pauli matrices represented by x1z1 and x2z2 are multiplied together + + :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: + 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. + """ + 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] + 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} as a submethod to applying CNOT to the tableau + + :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) + + def _apply_hadamard(self, instruction): + """ + Apply a hadamard gate on qubit defined in instruction + + :param instruction: pyquil abstract instruction. + """ + 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]] + + def _apply_phase(self, instruction): + """ + Apply the phase gate instruction ot the tableau + + :param instruction: pyquil abstract instruction. + """ + 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] + + 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) + + # 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 `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 + 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 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: # 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}) + # 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, :] 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 new file mode 100644 index 0000000..750a365 --- /dev/null +++ b/referenceqvm/stabilizer_utils.py @@ -0,0 +1,226 @@ +""" +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 +""" +from functools import reduce +from scipy.sparse import csc_matrix +import numpy as np +from pyquil.paulis import PauliTerm, sX, sZ, sY + + +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: new classical state and the coefficient it picked up. + """ + 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") + + # 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 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") + + 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 = [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) + + 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, + 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 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} + :return: state projected by stabilizers + """ + if num_qubits is None: + num_qubits = len(stabilizer_list) + + 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) + + for generator in stabilizer_list: + # (I + G(i)) / 2 + state += state_family_generator(state, generator) + state /= 2 + + normalization = (state.conj().T.dot(state)).toarray() + state /= np.sqrt(float(normalization.real)) + + 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_stabilizer_qvm.py b/referenceqvm/tests/test_stabilizer_qvm.py new file mode 100644 index 0000000..fbece46 --- /dev/null +++ b/referenceqvm/tests/test_stabilizer_qvm.py @@ -0,0 +1,400 @@ +""" +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, I +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 + + +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_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))])) + 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) + 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) + 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 + + +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) + + 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) + qvmstab._apply_phase(prog.instructions[0]) + true_stab = np.array([[1, 1, 0], + [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(): + """ + 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] + + # 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(): + """ + 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)).measure(0, 0) + qvmstab = QVM_Stabilizer() + 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) + 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) diff --git a/referenceqvm/tests/test_stabilizer_utils.py b/referenceqvm/tests/test_stabilizer_utils.py new file mode 100644 index 0000000..039aebe --- /dev/null +++ b/referenceqvm/tests/test_stabilizer_utils.py @@ -0,0 +1,294 @@ +""" +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.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 + + +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_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 + """ + 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, np.array(stabilizer_state.todense())) diff --git a/referenceqvm/unitary_generator.py b/referenceqvm/unitary_generator.py index 1cb5104..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,7 +417,7 @@ 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: + raise TypeError("Object not recognizable")