diff --git a/referenceqvm/api.py b/referenceqvm/api.py index abfeef9..422770c 100644 --- a/referenceqvm/api.py +++ b/referenceqvm/api.py @@ -4,10 +4,11 @@ from referenceqvm.gates import gate_matrix from referenceqvm.qvm_wavefunction import QVM_Wavefunction from referenceqvm.qvm_unitary import QVM_Unitary +from referenceqvm.qvm_density import QVM_Density def QVMConnection(type_trans='wavefunction', - gate_set=gate_matrix): + gate_set=gate_matrix, noise_model=None): """ Initialize a qvm of a particular type. The type corresponds to the type of transition the QVM can perform. @@ -31,6 +32,9 @@ def QVMConnection(type_trans='wavefunction', elif type_trans == 'unitary': qvm = QVM_Unitary(gate_set=gate_set) + elif type_trans == 'density': + qvm = QVM_Density(gate_set=gate_set, noise_model=noise_model) + else: raise TypeError("{} is not a valid QVM type.".format(type_trans)) diff --git a/referenceqvm/gates.py b/referenceqvm/gates.py index cdee908..6cf84be 100644 --- a/referenceqvm/gates.py +++ b/referenceqvm/gates.py @@ -226,3 +226,66 @@ def BARENCO(alpha, phi, theta): 'CZ': CZ } + +# noisy gates +def relaxation(p): + """ + Return the amplitude damping Kraus operators + """ + k0 = np.array([[1.0, 0.0], [0.0, np.sqrt(1 - p)]]) + k1 = np.array([[0.0, np.sqrt(p)], [0.0, 0.0]]) + return [k0, k1] + + +def dephasing(p): + """ + Return the phase damping Kraus operators + """ + k0 = np.eye(2)*np.sqrt(1 - p/2) + k1 = np.sqrt(p/2) * Z + return [k0, k1] + + +def depolarizing(p): + """ + Return the phase damping Kraus operators + """ + k0 = np.sqrt(1.0 - p)*I + k1 = np.sqrt(p/3.0)*X + k2 = np.sqrt(p/3.0)*Y + k3 = np.sqrt(p/3.0)*Z + return [k0, k1, k2, k3] + + +def phase_flip(p): + """ + Return the phase flip kraus operators + """ + k0 = np.sqrt(1 - p)*I + k1 = np.sqrt(p)*Z + return [k0, k1] + + +def bit_flip(p): + """ + Return the phase flip kraus operators + """ + k0 = np.sqrt(1 - p)*I + k1 = np.sqrt(p)*X + return [k0, k1] + + +def bitphase_flip(p): + """ + Return the phase flip kraus operators + """ + k0 = np.sqrt(1 - p)*I + k1 = np.sqrt(p)*Y + return [k0, k1] + +noise_gates = {'relaxation': relaxation, + 'dephasing': dephasing, + 'depolarizing': depolarizing, + 'phase_flip': phase_flip, + 'bit_flip': bit_flip, + 'bitphase_flip': bitphase_flip} \ No newline at end of file diff --git a/referenceqvm/qam.py b/referenceqvm/qam.py index 0e06ea2..c9e3507 100644 --- a/referenceqvm/qam.py +++ b/referenceqvm/qam.py @@ -89,6 +89,7 @@ def load_program(self, pyquil_program): 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") diff --git a/referenceqvm/qvm_density.py b/referenceqvm/qvm_density.py new file mode 100644 index 0000000..3cb9d96 --- /dev/null +++ b/referenceqvm/qvm_density.py @@ -0,0 +1,454 @@ +#!/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 Gates, and returns the +unitary resulting from the program evolution. +""" +import sys +import numpy as np +import scipy.sparse as sps + +from pyquil.quil import Program +from pyquil.quilbase import (Measurement, + UnaryClassicalInstruction, + BinaryClassicalInstruction, + ClassicalTrue, + ClassicalFalse, + ClassicalOr, + ClassicalNot, + ClassicalAnd, + ClassicalExchange, + ClassicalMove, + Jump, + JumpTarget, + JumpConditional, + JumpWhen, + JumpUnless, + Halt, Gate) + +from referenceqvm.qam import QAM +from referenceqvm.unitary_generator import lifted_gate, tensor_gates, value_get +from referenceqvm.gates import utility_gates, noise_gates + +INFINITY = float("inf") +"Used for infinite coherence times." + + +def sparse_trace(sparse_matrix): + dim = sparse_matrix.shape[0] + diagonal = [sparse_matrix[i, i] for i in range(dim)] + return np.sum(diagonal) + + +class NoiseModel(object): + + def __init__(self, T1=30e-6, T2=30e-6, gate_time_1q=50e-9, + gate_time_2q=150e-09, ro_fidelity=0.95, + depolarizing=None, bitflip=None, phaseflip=None, + bitphaseflip=None): + """ + Generate a noise model + + :param T1: Relaxation time. + :param T2: Dephasing time. + :param gate_time_1q: 1-qubit gate time. + :param gate_time_2q: 2-qubit gate time. + :param ro_fidelity: read-out fidelity (constructs a symmetric readout + confusion matrix> + """ + self.T1 = T1 + self.T2 = T2 + self.gate_time_1q = gate_time_1q + self.gate_time_2q = gate_time_2q + self.gate_time_3q = float(200E-9) + self.gate_time_4q = float(200E-9) + self.gate_time_5q = float(200E-9) + self.ro_fidelity = ro_fidelity + + # alternative error models + self.depolarizing = depolarizing + self.bitflip = bitflip + self.phaseflip = phaseflip + self.bitphaseflip = bitphaseflip + + +class QVM_Density(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, + density=None, noise_model=None): + """ + Subclassed from QAM this is a pure QVM. + """ + super(QVM_Density, 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._density = density + self.noise_model = noise_model + self.all_inst = True + + def measurement(self, qubit_index): + """ + Perform measurement + + :param qubit_index: qubit to measure + :return: + """ + measure_0 = lifted_gate(qubit_index, utility_gates['P0'], + self.num_qubits) + prob_zero = sparse_trace(measure_0.dot(self._density)) + # generate random number to 'roll' for measurement + if np.random.random() < prob_zero: + # decohere state using the measure_0 operator + unitary = measure_0.dot( + sps.eye(2 ** self.num_qubits) / np.sqrt(prob_zero)) + measured_val = 0 + else: # measure one + measure_1 = lifted_gate(qubit_index, utility_gates['P1'], + self.num_qubits) + unitary = measure_1.dot( + sps.eye(2 ** self.num_qubits) / np.sqrt(1 - prob_zero)) + measured_val = 1 + + return measured_val, unitary + + def _transition(self, instruction): + """ + Implements a transition on the density-qvm. + + :param Gate instruction: QuilAction gate to be implemented + """ + if isinstance(instruction, Measurement): + # perform measurement and modify wf in-place + t_qbit = value_get(instruction.qubit) + t_cbit = value_get(instruction.classical_reg) + measured_val, unitary = self.measurement(t_qbit) + self._density = unitary.dot(self._density).dot(np.conj(unitary.T)) + + # load measured value into classical bit destination + self.classical_memory[t_cbit] = measured_val + self.program_counter += 1 + + elif isinstance(instruction, Gate): + # apply Gate or DefGate + unitary = tensor_gates(self.gate_set, self.defgate_set, instruction, self.num_qubits) + self._density = unitary.dot(self._density).dot(np.conj(unitary).T) + 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 + + elif isinstance(instruction, JumpConditional): + # JumpConditional; check classical reg + cond = self.classical_memory[value_get(instruction.condition)] + dest_index = self.find_label(instruction.target) + if isinstance(instruction, JumpWhen): + jump_if_cond = True + elif isinstance(instruction, JumpUnless): + jump_if_cond = False + else: + raise TypeError("Invalid JumpConditional") + + if not (cond ^ jump_if_cond): + # jumping: set prog counter to JumpTarget + self.program_counter = dest_index + else: + # not jumping: hop over this JumpConditional + self.program_counter += 1 + + elif isinstance(instruction, UnaryClassicalInstruction): + # UnaryClassicalInstruction; set classical reg + target_ind = value_get(instruction.target) + old = self.classical_memory[target_ind] + if isinstance(instruction, ClassicalTrue): + new = True + elif isinstance(instruction, ClassicalFalse): + new = False + elif isinstance(instruction, ClassicalNot): + new = not old + else: + raise TypeError("Invalid UnaryClassicalInstruction") + + self.classical_memory[target_ind] = new + self.program_counter += 1 + + elif isinstance(instruction, BinaryClassicalInstruction): + # BinaryClassicalInstruction; compute and set classical reg + left_ind = value_get(instruction.left) + left_val = self.classical_memory[left_ind] + right_ind = value_get(instruction.right) + right_val = self.classical_memory[right_ind] + if isinstance(instruction, ClassicalAnd): + # compute LEFT AND RIGHT, set RIGHT to the result + self.classical_memory[right_ind] = left_val & right_val + elif isinstance(instruction, ClassicalOr): + # compute LEFT OR RIGHT, set RIGHT to the result + self.classical_memory[right_ind] = left_val | right_val + elif isinstance(instruction, ClassicalMove): + # move LEFT to RIGHT + self.classical_memory[right_ind] = left_val + elif isinstance(instruction, ClassicalExchange): + # exchange LEFT and RIGHT + self.classical_memory[left_ind] = right_val + self.classical_memory[right_ind] = left_val + else: + raise TypeError("Invalid BinaryClassicalInstruction") + + self.program_counter += 1 + elif isinstance(instruction, Halt): + # do nothing; set program_counter to end of program + self.program_counter = len(self.program) + else: + raise TypeError("Invalid / unsupported instruction type: {}\n" + "Currently supported: unary/binary classical " + "instructions, control flow (if/while/jumps), " + "measurements, and gates/defgates.".format(type(instruction))) + + def _pre(self, instruction): + """ + Pre-hook for state-machine execution + + Useful for applying error prior to measurement + + :param instruction: Instruction to apply to the state of the + state-machine + :return: None + """ + if self.noise_model is None: + return + + else: + + if isinstance(instruction, Measurement): + # perform bitflip with probability associated with ro_fidelity + # 1 - is because we cite fidelity not infidelity + ro_prob = 1 - self.noise_model.ro_fidelity + bitflip_kraus_ops = noise_gates['bitphase_flip'](ro_prob) + for iqubit in range(self.num_qubits): + for kop in bitflip_kraus_ops: + kraus_op = lifted_gate(iqubit, kop, self.num_qubits) + self._density += kraus_op.dot(self._density) + + def _post(self, instruction): + """ + Post-hook for state-machine execution + + useful for applying error after the gate + + :param instruction: Instruction to apply to the state of the + state-machine + :return: None + """ + if self.noise_model is None: + return + else: + if self.noise_model.T1 != INFINITY: + if isinstance(instruction, Gate): + # perform kraus operator for relaxation + n_qubit_span = len(instruction.get_qubits()) + if n_qubit_span == 1: + p = 1.0 - np.exp(-self.noise_model.gate_time_1q/self.noise_model.T1) + elif n_qubit_span == 2: + p = 1.0 - np.exp(-self.noise_model.gate_time_2q/self.noise_model.T1) + else: + # optional to expand to larger number of qubits + p = 1.0 - np.exp(-self.noise_model.gate_time_2q/self.noise_model.T1) + + if not np.isclose(p, 1.0): + relaxation_ops = noise_gates['relaxation'](p) + + for i in range(self.num_qubits): + self._density = self.apply_kraus(relaxation_ops, i) + + if self.noise_model.T2 != INFINITY: + # perform kraus operator for relaxation + n_qubit_span = len(instruction.get_qubits()) + if n_qubit_span == 1: + p = 1.0 - np.exp(-self.noise_model.gate_time_1q/self.noise_model.T2) + elif n_qubit_span == 2: + p = 1.0 - np.exp(-self.noise_model.gate_time_2q/self.noise_model.T2) + else: + # optional to expand to larger number of qubits + p = 1.0 - np.exp(-self.noise_model.gate_time_2q/self.noise_model.T2) + + if not np.isclose(p, 1.0): + dephasing_ops = noise_gates['dephasing'](p) + + for i in range(self.num_qubits): + self._density = self.apply_kraus(dephasing_ops, i) + + if self.noise_model.depolarizing is not None: + if not isinstance(self.noise_model.depolarizing, float): + raise TypeError("depolarizing noise value should be None or a float") + + depolarizing_kraus_ops = noise_gates['depolarizing'](self.noise_model.depolarizing) + + for i in range(self.num_qubits): + self._density = self.apply_kraus(depolarizing_kraus_ops, i) + + if self.noise_model.bitflip is not None: + if not isinstance(self.noise_model.bitflip, float): + raise TypeError("bitflip noise value should be None or a float") + + bitflip_kraus_ops = noise_gates['bit_flip'](self.noise_model.bitflip) + for i in range(self.num_qubits): + self._density = self.apply_kraus(bitflip_kraus_ops, i) + + if self.noise_model.phaseflip is not None: + if not isinstance(self.noise_model.phaseflip, float): + raise TypeError("phaseflip noise value should be None or a float") + + phaseflip_kraus_ops = noise_gates['phase_flip'](self.noise_model.phaseflip) + for i in range(self.num_qubits): + self._density = self.apply_kraus(phaseflip_kraus_ops, i) + + if self.noise_model.bitphaseflip is not None: + if not isinstance(self.noise_model.phaseflip, float): + raise TypeError("bitphaseflip noise value should be None or a float") + + bitphaseflip_kraus_ops = noise_gates['bitphase_flip'](self.noise_model.bitphaseflip) + for i in range(self.num_qubits): + self._density = self.apply_kraus(bitphaseflip_kraus_ops, i) + + return + + def apply_kraus(self, operators, index): + """ + Apply Kraus operators to density + """ + output_density = sps.csc_matrix((2**self.num_qubits, + 2**self.num_qubits), + dtype=complex) + + for op in operators: + large_op = sps.csc_matrix(lifted_gate(index, op, self.num_qubits)) + output_density += large_op.dot(self._density).dot(np.conj(large_op.T)) + return output_density + + def transition(self, instruction): + """ + Implements a transition on the density-qvm with pre and post hooks + + :param Gate instruction: QuilAction gate to be implemented + """ + self._pre(instruction) + self._transition(instruction) + self._post(instruction) + + def density(self, pyquil_program): + """ + Return the density matrix 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 density matrix corresponding to the output of the program. + :rtype: np.array + """ + + # load program + self.load_program(pyquil_program) + # setup density + N = 2 ** self.num_qubits + self._density = sps.csc_matrix(([1.0], ([0], [0])), shape=(N, N)) + # evolve unitary + self.kernel() + return self._density + + 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 + """ + raise NotImplementedError() + + def run(self, pyquil_program, classical_addresses=None, trials=1): + """ + Run a pyQuil program multiple times, accumulating the values deposited + in a list of classical addresses. + + This uses numpy's inverse sampling method to calculate bit string + outcomes + + :param Program pyquil_program: A pyQuil program. + :param list classical_addresses: A list of classical addresses. + This is ignored but kept to have + similar input as Forest QVM. + :param int trials: Number of shots to collect. + + :return: A list of lists of bits. Each sublist corresponds to the + values in `classical_addresses`. + :rtype: list + """ + results = [] + for trial in range(trials): + _ = self.density(pyquil_program) + results.append(self.classical_memory[classical_addresses]) + return list(results) + + def run_and_measure(self, pyquil_program, trials=1): + """ + Run and measure. converts to run + + :param pyquil_program: + :param trials: + :return: + """ + density = self.density(pyquil_program) + probabilities = [density[i, i].real for i in range(2 ** self.num_qubits)] + + results_as_integers = np.random.choice(2 ** self.num_qubits, + size=trials, + p=probabilities) + + results = list(map(lambda x: + list(map(int, np.binary_repr(x, width=self.num_qubits))), + results_as_integers)) + return results + + diff --git a/referenceqvm/tests/test_density.py b/referenceqvm/tests/test_density.py new file mode 100644 index 0000000..c42ec49 --- /dev/null +++ b/referenceqvm/tests/test_density.py @@ -0,0 +1,377 @@ +import numpy as np +import scipy.sparse as sps +from math import pi +from pyquil.quil import Program +from pyquil.gates import H as Hgate +from pyquil.gates import CNOT as CNOTgate +from pyquil.gates import CZ as CZgate +from pyquil.gates import X as Xgate +from pyquil.gates import RX as RXgate +from pyquil.gates import RY as RYgate +from pyquil.gates import RZ as RZgate +from pyquil.gates import I as Igate +from pyquil.gates import PHASE as PHASEgate +from referenceqvm.api import QVMConnection +from referenceqvm.qvm_density import QVM_Density, NoiseModel, INFINITY, lifted_gate +from referenceqvm.gates import noise_gates, X, Y, Z, I, P0, P1, RY + + +def random_1q_density(): + state = np.random.random(2) + 1j*np.random.random() + normalization = np.conj(state).T.dot(state) + state /= np.sqrt(normalization) + state = state.reshape((-1, 1)) + + rho = state.dot(np.conj(state).T) + assert np.isclose(np.trace(rho), 1.0) + assert np.allclose(rho, np.conj(rho).T) + return rho + + +def test_initialize(): + qvm = QVMConnection(type_trans='density') + assert isinstance(qvm, QVM_Density) + + +def test_qaoa_density(): + wf_true = [0.00167784 + 1.00210180e-05*1j, 0.50000000 - 4.99997185e-01*1j, + 0.50000000 - 4.99997185e-01*1j, 0.00167784 + 1.00210180e-05*1j] + wf_true = np.reshape(np.array(wf_true), (4, 1)) + rho_true = np.dot(wf_true, np.conj(wf_true).T) + prog = Program() + prog.inst([RYgate(np.pi/2)(0), RXgate(np.pi)(0), + RYgate(np.pi/2)(1), RXgate(np.pi)(1), + CNOTgate(0, 1), RXgate(-np.pi/2)(1), RYgate(4.71572463191)(1), + RXgate(np.pi/2)(1), CNOTgate(0, 1), + RXgate(-2*2.74973750579)(0), RXgate(-2*2.74973750579)(1)]) + + qvm = QVMConnection(type_trans='density') + wf_rho = qvm.density(prog) + assert np.isclose(rho_true, wf_rho.todense()).all() + + +def test_larger_qaoa_density(): + square_qaoa_circuit = [Hgate(0), Hgate(1), Hgate(2), Hgate(3), + Xgate(0), + PHASEgate(0.3928244130249029)(0), + Xgate(0), + PHASEgate(0.3928244130249029)(0), + CNOTgate(0, 1), + RZgate(0.78564882604980579)(1), + CNOTgate(0, 1), + Xgate(0), + PHASEgate(0.3928244130249029)(0), + Xgate(0), + PHASEgate(0.3928244130249029)(0), + CNOTgate(0, 3), + RZgate(0.78564882604980579)(3), + CNOTgate(0, 3), + Xgate(0), + PHASEgate(0.3928244130249029)(0), + Xgate(0), + PHASEgate(0.3928244130249029)(0), + CNOTgate(1, 2), + RZgate(0.78564882604980579)(2), + CNOTgate(1, 2), + Xgate(0), + PHASEgate(0.3928244130249029)(0), + Xgate(0), + PHASEgate(0.3928244130249029)(0), + CNOTgate(2, 3), + RZgate(0.78564882604980579)(3), + CNOTgate(2, 3), + Hgate(0), + RZgate(-0.77868204192240842)(0), + Hgate(0), + Hgate(1), + RZgate(-0.77868204192240842)(1), + Hgate(1), + Hgate(2), + RZgate(-0.77868204192240842)(2), + Hgate(2), + Hgate(3), + RZgate(-0.77868204192240842)(3), + Hgate(3)] + + prog = Program(square_qaoa_circuit) + qvm = QVMConnection(type_trans='density') + rho_test = qvm.density(prog) + wf_true = np.array([8.43771693e-05-0.1233845*1j, -1.24927731e-01+0.00329533*1j, + -1.24927731e-01+0.00329533*1j, + -2.50040954e-01+0.12661547*1j, + -1.24927731e-01+0.00329533*1j, -4.99915497e-01-0.12363516*1j, + -2.50040954e-01+0.12661547*1j, -1.24927731e-01+0.00329533*1j, + -1.24927731e-01+0.00329533*1j, -2.50040954e-01+0.12661547*1j, + -4.99915497e-01-0.12363516*1j, -1.24927731e-01+0.00329533*1j, + -2.50040954e-01+0.12661547*1j, -1.24927731e-01+0.00329533*1j, + -1.24927731e-01+0.00329533*1j, + 8.43771693e-05-0.1233845*1j]) + + wf_true = np.reshape(wf_true, (2**4, 1)) + rho_true = np.dot(wf_true, np.conj(wf_true).T) + assert np.isclose(rho_test.todense(), rho_true).all() + + +def get_compiled_prog(theta): + return Program([ + RZgate(-pi/2, 0), + RXgate(-pi/2, 0), + RZgate(-pi/2, 1), + RXgate( pi/2, 1), + CZgate(1, 0), + RZgate(-pi/2, 1), + RXgate(-pi/2, 1), + RZgate(theta, 1), + RXgate( pi/2, 1), + CZgate(1, 0), + RXgate( pi/2, 0), + RZgate( pi/2, 0), + RZgate(-pi/2, 1), + RXgate( pi/2, 1), + RZgate(-pi/2, 1), + ]) + + +def test_kraus_t1_normalization(): + kraus_ops = noise_gates['relaxation'](0.75) + total = np.zeros((2, 2), dtype=complex) + for kop in kraus_ops: + total += np.conj(kop.T).dot(kop) + assert np.allclose(total, np.eye(2)) + + +def test_kraus_t2_normalization(): + kraus_ops = noise_gates['dephasing'](0.75) + total = np.zeros((2, 2), dtype=complex) + for kop in kraus_ops: + total += np.conj(kop.T).dot(kop) + assert np.allclose(total, np.eye(2)) + + +def test_kraus_depolarizing_normalization(): + kraus_ops = noise_gates['depolarizing'](0.75) + total = np.zeros((2, 2), dtype=complex) + for kop in kraus_ops: + total += np.conj(kop.T).dot(kop) + assert np.allclose(total, np.eye(2)) + + +def test_kraus_bitflip_normalization(): + kraus_ops = noise_gates['bit_flip'](0.75) + total = np.zeros((2, 2), dtype=complex) + for kop in kraus_ops: + total += np.conj(kop.T).dot(kop) + assert np.allclose(total, np.eye(2)) + + +def test_kraus_phaseflip_normalization(): + kraus_ops = noise_gates['phase_flip'](0.75) + total = np.zeros((2, 2), dtype=complex) + for kop in kraus_ops: + total += np.conj(kop.T).dot(kop) + assert np.allclose(total, np.eye(2)) + + +def test_kraus_bitphaseflip_normalization(): + kraus_ops = noise_gates['bitphase_flip'](0.75) + total = np.zeros((2, 2), dtype=complex) + for kop in kraus_ops: + total += np.conj(kop.T).dot(kop) + assert np.allclose(total, np.eye(2)) + + +def test_kraus_application_bitflip(): + qvm = QVMConnection(type_trans='density') + initial_density = random_1q_density() + qvm._density = sps.csc_matrix(initial_density) + qvm.num_qubits = 1 + p = 0.372 + kraus_ops = noise_gates['bit_flip'](p) + final_density = (1 - p) * initial_density + p * X.dot(initial_density).dot(X) + test_density = qvm.apply_kraus(kraus_ops, 0) + assert np.allclose(test_density.todense(), final_density) + + +def test_kraus_application_phaseflip(): + qvm = QVMConnection(type_trans='density') + initial_density = random_1q_density() + qvm._density = sps.csc_matrix(initial_density) + qvm.num_qubits = 1 + p = 0.372 + kraus_ops = noise_gates['phase_flip'](p) + final_density = (1 - p) * initial_density + p * Z.dot(initial_density).dot(Z) + test_density = qvm.apply_kraus(kraus_ops, 0) + assert np.allclose(test_density.todense(), final_density) + + +def test_kraus_application_bitphaseflip(): + qvm = QVMConnection(type_trans='density') + initial_density = random_1q_density() + qvm._density = sps.csc_matrix(initial_density) + qvm.num_qubits = 1 + p = 0.372 + kraus_ops = noise_gates['bitphase_flip'](p) + final_density = (1 - p) * initial_density + p * Y.dot(initial_density).dot(Y) + test_density = qvm.apply_kraus(kraus_ops, 0) + assert np.allclose(test_density.todense(), final_density) + + +def test_kraus_application_relaxation(): + qvm = QVMConnection(type_trans='density') + rho = random_1q_density() + qvm._density = sps.csc_matrix(rho) + qvm.num_qubits = 1 + p = 0.372 + kraus_ops = noise_gates['relaxation'](p) + final_density = np.array([[rho[0, 0] + rho[1, 1] * p, np.sqrt(1 - p) * rho[0, 1]], + [np.sqrt(1 - p) * rho[1, 0], (1 - p) * rho[1, 1]]]) + test_density = qvm.apply_kraus(kraus_ops, 0) + assert np.allclose(test_density.todense(), final_density) + + +def test_kraus_application_dephasing(): + qvm = QVMConnection(type_trans='density') + rho = random_1q_density() + qvm._density = sps.csc_matrix(rho) + qvm.num_qubits = 1 + p = 0.372 + kraus_ops = noise_gates['dephasing'](p) + final_density = np.array([[rho[0, 0], (1 - p) * rho[0, 1]], + [(1 - p) * rho[1, 0], rho[1, 1]]]) + test_density = qvm.apply_kraus(kraus_ops, 0) + assert np.allclose(test_density.todense(), final_density) + + +def test_kraus_application_depolarizing(): + qvm = QVMConnection(type_trans='density') + rho = random_1q_density() + qvm._density = sps.csc_matrix(rho) + qvm.num_qubits = 1 + p = 0.372 + kraus_ops = noise_gates['depolarizing'](p) + final_density = (1 - p) * rho + (p / 3) * (X.dot(rho).dot(X) + + Y.dot(rho).dot(Y) + + Z.dot(rho).dot(Z)) + test_density = qvm.apply_kraus(kraus_ops, 0) + assert np.allclose(test_density.todense(), final_density) + + +def test_kraus_compound_T1T2_application(): + qvm = QVMConnection(type_trans='density') + rho = random_1q_density() + qvm._density = sps.csc_matrix(rho) + qvm.num_qubits = 1 + p1 = 0.372 + p2 = 0.45 + kraus_ops_1 = noise_gates['relaxation'](p1) + kraus_ops_2 = noise_gates['dephasing'](p2) + final_density = np.array([[rho[0, 0] + rho[1, 1] * p1, (1 - p2 ) * np.sqrt(1 - p1) * rho[0, 1]], + [(1 - p2) * np.sqrt(1 - p1) * rho[1, 0], (1 - p1) * rho[1, 1]]]) + qvm._density = qvm.apply_kraus(kraus_ops_1, 0) + qvm._density = qvm.apply_kraus(kraus_ops_2, 0) + assert np.allclose(qvm._density.todense(), final_density) + + +def test_kraus_through_qvm_t1(): + noise = NoiseModel(T2=INFINITY, ro_fidelity=1.0) + qvm = QVMConnection(type_trans='density', noise_model=noise) + rho = random_1q_density() + identity_program = Program().inst(Igate(0)) + qvm._density = sps.csc_matrix(rho) + qvm.num_qubits = 1 + qvm.load_program(identity_program) + qvm.kernel() + p = 1 - np.exp(-noise.gate_time_1q/noise.T1) + final_density = np.array([[rho[0, 0] + rho[1, 1] * p, np.sqrt(1 - p) * rho[0, 1]], + [np.sqrt(1 - p) * rho[1, 0], (1 - p) * rho[1, 1]]]) + assert np.allclose(final_density, qvm._density.todense()) + + +def test_kraus_through_qvm_t2(): + noise = NoiseModel(T1=INFINITY, ro_fidelity=1.0) + qvm = QVMConnection(type_trans='density', noise_model=noise) + rho = random_1q_density() + identity_program = Program().inst(Igate(0)) + qvm._density = sps.csc_matrix(rho) + qvm.num_qubits = 1 + qvm.load_program(identity_program) + qvm.kernel() + p = 1 - np.exp(-noise.gate_time_1q/noise.T2) + final_density = np.array([[rho[0, 0], (1 - p) * rho[0, 1]], + [(1 - p) * rho[1, 0], rho[1, 1]]]) + assert np.allclose(final_density, qvm._density.todense()) + + +def test_kraus_through_qvm_t1t2(): + noise = NoiseModel(ro_fidelity=1.0) + qvm = QVMConnection(type_trans='density', noise_model=noise) + rho = random_1q_density() + identity_program = Program().inst(Igate(0)) + qvm._density = sps.csc_matrix(rho) + qvm.num_qubits = 1 + qvm.load_program(identity_program) + qvm.kernel() + p = 1 - np.exp(-noise.gate_time_1q/noise.T2) + p1 = 1 - np.exp(-noise.gate_time_1q/noise.T1) + final_density = np.array([[rho[0, 0] + p1 * rho[1, 1], np.sqrt(1 - p1) * (1 - p) * rho[0, 1]], + [np.sqrt(1 - p1) * (1 - p) * rho[1, 0], (1 - p1) * rho[1, 1]]]) + assert np.allclose(final_density, qvm._density.todense()) + + +def multi_qubit_decay_bellstate(): + """ + Test multiqubit decay + """ + program = Program().inst([RYgate(np.pi/3)(0), CNOTgate(0, 1)]) + noise = NoiseModel(ro_fidelity=1.0) + qvm = QVMConnection(type_trans='density', noise_model=noise) + + initial_density = np.zeros((4, 4), dtype=complex) + initial_density[0, 0] = 1.0 + cnot_01 = np.kron(I, P0) + np.kron(X, P1) + + p1 = 1 - np.exp(-noise.gate_time_1q/noise.T1) + p2 = 1 - np.exp(-noise.gate_time_1q/noise.T2) + + kraus_ops_1 = noise_gates['relaxation'](p1) + kraus_ops_2 = noise_gates['dephasing'](p2) + + gate_1 = np.kron(np.eye(2), RY(np.pi/3)) + + state = gate_1.dot(initial_density).dot(np.conj(gate_1).T) + for ii in range(2): + new_density = np.zeros_like(state) + for kop in kraus_ops_1: + operator = lifted_gate(ii, kop, 2).todense() + new_density += operator.dot(state).dot(np.conj(operator).T) + state = new_density + + for ii in range(2): + new_density = np.zeros_like(state) + for kop in kraus_ops_2: + operator = lifted_gate(ii, kop, 2).todense() + new_density += operator.dot(state).dot(np.conj(operator).T) + state = new_density + + state = cnot_01.dot(state).dot(cnot_01.T) + p1 = 1 - np.exp(-noise.gate_time_2q/noise.T1) + p2 = 1 - np.exp(-noise.gate_time_2q/noise.T2) + kraus_ops_1 = noise_gates['relaxation'](p1) + kraus_ops_2 = noise_gates['dephasing'](p2) + + for ii in range(2): + new_density = np.zeros_like(state) + for kop in kraus_ops_1: + operator = lifted_gate(ii, kop, 2).todense() + new_density += operator.dot(state).dot(np.conj(operator).T) + state = new_density + + for ii in range(2): + new_density = np.zeros_like(state) + for kop in kraus_ops_2: + operator = lifted_gate(ii, kop, 2).todense() + new_density += operator.dot(state).dot(np.conj(operator).T) + state = new_density + + density = qvm.density(program) + assert np.allclose(density.todense(), state) diff --git a/referenceqvm/tests/test_measurement.py b/referenceqvm/tests/test_measurement.py index 39084b9..7574574 100644 --- a/referenceqvm/tests/test_measurement.py +++ b/referenceqvm/tests/test_measurement.py @@ -32,3 +32,5 @@ def test_biased_coin(qvm): results = qvm.run_and_measure(prog, qubits=[0], trials=samples) coin_bias = sum(map(lambda x: x[0], results)) / float(samples) assert np.isclose(coin_bias, 0.25, atol=0.05, rtol=0.05) + + diff --git a/referenceqvm/tests/test_sampling.py b/referenceqvm/tests/test_sampling.py new file mode 100644 index 0000000..6f39e2c --- /dev/null +++ b/referenceqvm/tests/test_sampling.py @@ -0,0 +1,54 @@ +""" +Testing sampling of a density matrix +""" +import numpy as np +from itertools import product + +from referenceqvm.api import QVMConnection +from pyquil.quil import Program +from pyquil.gates import H, CNOT, RX + + +def test_sample_coin(): + # Flip a coin with a density matrix. + prog = Program().inst(H(0)) + qvm = QVMConnection(type_trans='density') + + samples = 10000 + results = qvm.run_and_measure(prog, trials=samples) + coin_bias = sum(map(lambda x: x[0], results))/float(samples) + assert np.isclose(coin_bias, 0.5, atol=0.05, rtol=0.05) + + +def test_sample_bell(): + # Sample a bell state + prog = Program().inst(H(0), CNOT(0, 1)) + qvm = QVMConnection(type_trans='density') + + samples = 100000 + results = qvm.run_and_measure(prog, trials=samples) + for rr in results: + assert rr[0] == rr[1] + + bias_pair = sum(map(lambda x: x[0], results))/float(samples) + assert np.isclose(bias_pair, 0.5, atol=0.05, rtol=0.05) + + +def test_biased_coin(): + # sample from a %75 head and 25% tails coin + prog = Program().inst(RX(np.pi/3)(0)) + qvm = QVMConnection(type_trans='density') + samples = 100000 + results = qvm.run_and_measure(prog, trials=samples) + coin_bias = sum(map(lambda x: x[0], results))/float(samples) + assert np.isclose(coin_bias, 0.25, atol=0.05, rtol=0.05) + + +def test_measurement(): + qvm = QVMConnection(type_trans='density') + prog = Program().inst(H(0)).measure(0, [0]).measure(0, [1]) + samples = 10000 + result = qvm.run(prog, [0, 1], trials=samples) + assert all(map(lambda x: x[0] == x[1], result)) + bias = sum(map(lambda x: x[0], result))/float(samples) + assert np.isclose(0.5, bias, atol=0.05, rtol=0.05) diff --git a/referenceqvm/unitary_generator.py b/referenceqvm/unitary_generator.py index 78af8ab..1cb5104 100644 --- a/referenceqvm/unitary_generator.py +++ b/referenceqvm/unitary_generator.py @@ -21,6 +21,7 @@ Note: uses SciPy sparse diagonal (DIA) representation to increase space and timeefficiency. """ +import numpy as np from collections import Sequence from numbers import Integral @@ -355,7 +356,6 @@ def tensor_gates(gate_set, defgate_set, pyquil_gate, num_qubits): gate = apply_gate(dict_check[pyquil_gate.name], args, num_qubits) - return gate