From 80658dfe0548ab3f4bd9b1e85d3ce99aa20354bb Mon Sep 17 00:00:00 2001 From: Boxi Li Date: Sun, 22 Aug 2021 17:33:28 +0200 Subject: [PATCH 1/2] add truncation and expansion with qubit state --- src/qutip_qip/qubits.py | 144 +++++++++++++++++++++++++++++++++++++++- tests/test_qubits.py | 83 +++++++++++++++++++++-- 2 files changed, 220 insertions(+), 7 deletions(-) diff --git a/src/qutip_qip/qubits.py b/src/qutip_qip/qubits.py index f14b7cb8b..70568f00c 100644 --- a/src/qutip_qip/qubits.py +++ b/src/qutip_qip/qubits.py @@ -31,13 +31,15 @@ # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ############################################################################### -__all__ = ['qubit_states'] +__all__ = [ + 'qubit_states', 'truncate_to_qubit_state', 'expand_qubit_state'] +import qutip from qutip import (tensor, basis) +import numpy as np from numpy import sqrt - def qubit_states(N=1, states=[0]): """ Function to define initial state of the qubits. @@ -64,3 +66,141 @@ def qubit_states(N=1, states=[0]): return tensor([alpha * basis(2, 1) + sqrt(1 - alpha**2) * basis(2, 0) for alpha in state_list]) + + +def _find_reduced_indices(dims): + """Find the forresponding indices for a given qu""" + if len(dims) == 1: + return np.array([0, 1], dtype=np.intp) + else: + rest_reduced_indices = _find_reduced_indices(dims[1:]) + rest_full_dims = np.product(dims[1:]) + # [indices + 0, indices + a] + reduced_indices = ( + np.array([[0], [rest_full_dims]], dtype=np.intp) + + np.array([rest_reduced_indices, rest_reduced_indices], + dtype=np.intp) + ) + return reduced_indices.reshape(reduced_indices.size) + + +def truncate_to_qubit_state(state): + """ + Truncate a given quantum state into the computational subspace, + i.e., each subspace is a 2-level system such as + ``dims=[[2, 2, 2...],[2, 2, 2...]]``. + The non-computational state will be discarded. + Notice that the trace truncated state is in general small than 1. + + Parameters + ---------- + state : :obj:`qutip.Qobj` + The input quantum state, either a ket, a bra or a square operator. + + Returns + ------- + truncated_state : :obj:`qutip.Qobj` + The truncated state. + + Examples + -------- + >>> import qutip + >>> from qutip_qip.qubits import truncate_to_qubit_state + >>> state = qutip.rand_ket(6, dims=[[2, 3], [1, 1]], seed=0) + >>> state # doctest: +NORMALIZE_WHITESPACE + Quantum object: dims = [[2, 3], [1, 1]], shape = (6, 1), type = ket + Qobj data = + [[-0.10369056+0.02570495j] + [ 0.34852171+0.06053252j] + [-0.05552249+0.37861125j] + [ 0.41247492-0.38160681j] + [ 0.25951892-0.36729024j] + [ 0.12978757-0.42681426j]] + >>> truncate_to_qubit_state(state) # doctest: +NORMALIZE_WHITESPACE + Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket + Qobj data = + [[-0.10369056+0.02570495j] + [ 0.34852171+0.06053252j] + [ 0.41247492-0.38160681j] + [ 0.25951892-0.36729024j]] + """ + + if state.isbra: + dims = state.dims[1] + reduced_dims = [[1] * len(dims), [2] * len(dims)] + elif state.isket: + dims = state.dims[0] + reduced_dims = [[2] * len(dims), [1] * len(dims)] + elif state.isoper and state.dims[0] == state.dims[1]: + dims = state.dims[0] + reduced_dims = [[2] * len(dims), [2] * len(dims)] + else: + raise ValueError("Can only truncate bra, ket or square operator.") + reduced_indices = _find_reduced_indices(dims) + # See NumPy Advanced Indexing. Choose the corresponding rows and columns. + zero_slice = np.array([0], dtype=np.intp) + reduced_indices1 = reduced_indices if not state.isbra else zero_slice + reduced_indices2 = reduced_indices if not state.isket else zero_slice + output = state[reduced_indices1[:, np.newaxis], reduced_indices2] + return qutip.Qobj(output, dims=reduced_dims) + + +def expand_qubit_state(state, dims): + """ + Expand a given qubit state into the a quantum state into a + multi-dimensional quantum state. + The additional matrix entries are filled with zeros. + + Parameters + ---------- + state : :obj:`qutip.Qobj` + The input quantum state, either a ket, a bra or a square operator. + + Returns + ------- + expanded_state : :obj:`qutip.Qobj` + The expanded state. + + Examples + -------- + >>> import qutip + >>> from qutip_qip.qubits import expand_qubit_state + >>> state = qutip.rand_ket(4, dims=[[2, 2], [1, 1]], seed=0) + >>> state # doctest: +NORMALIZE_WHITESPACE + Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket + Qobj data = + [[ 0.22362331-0.09566496j] + [-0.11702026+0.60050111j] + [ 0.15751346+0.71069216j] + [ 0.06879596-0.1786583j ]] + >>> expand_qubit_state(state, [3, 2]) # doctest: +NORMALIZE_WHITESPACE + Quantum object: dims = [[3, 2], [1, 1]], shape = (6, 1), type = ket + Qobj data = + [[ 0.22362331-0.09566496j] + [-0.11702026+0.60050111j] + [ 0.15751346+0.71069216j] + [ 0.06879596-0.1786583j ] + [ 0. +0.j ] + [ 0. +0.j ]] + """ + if state.isbra: + reduced_dims = state.dims[1] + full_dims = [[1] * len(dims), dims] + elif state.isket: + reduced_dims = state.dims[0] + full_dims = [dims, [1] * len(dims)] + elif state.isoper and state.dims[0] == state.dims[1]: + reduced_dims = state.dims[0] + full_dims = [dims, dims] + else: + raise ValueError("Can only expand bra, ket or square operator.") + if not all([d == 2 for d in reduced_dims]): + raise ValueError("The input state is not a qubit state.") + reduced_indices = _find_reduced_indices(dims) + + zero_slice = np.array([0], dtype=np.intp) + output = np.zeros([np.product(d) for d in full_dims], dtype=complex) + reduced_indices1 = reduced_indices if not state.isbra else zero_slice + reduced_indices2 = reduced_indices if not state.isket else zero_slice + output[reduced_indices1[:, np.newaxis], reduced_indices2] = state + return qutip.Qobj(output, dims=full_dims) diff --git a/tests/test_qubits.py b/tests/test_qubits.py index 757b764b9..5e6250dff 100644 --- a/tests/test_qubits.py +++ b/tests/test_qubits.py @@ -30,11 +30,13 @@ # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ############################################################################### - - from numpy.testing import assert_, run_module_suite -from qutip_qip.qubits import qubit_states +import numpy as np +import pytest +import qutip from qutip import (tensor, basis) +from qutip_qip.qubits import ( + qubit_states, truncate_to_qubit_state, expand_qubit_state) class TestQubits: @@ -57,6 +59,77 @@ def testQubitStates(self): psi01_b = qubit_states(N=2, states=[0, 1]) assert_(psi01_a == psi01_b) + @pytest.mark.parametrize( + "state, full_dims", + [(qutip.rand_dm(18, dims=[[3, 2, 3], [3, 2, 3]]), [3, 2, 3]), + (qutip.rand_ket(18, dims=[[2, 3, 3], [1, 1, 1]]), [2, 3, 3]), + (qutip.Qobj( + qutip.rand_ket(18).full().transpose(), + dims=[[1, 1, 1], [3, 2, 3]]), + [3, 2, 3] + ) + ] + ) + def test_state_truncation(self, state, full_dims): + reduced_state = truncate_to_qubit_state(state) + for _ in range(5): + ind1 = np.random.choice([0, 1], len(full_dims)) + ind2 = np.random.choice([0, 1], len(full_dims)) + ind_red1 = np.ravel_multi_index(ind1, [2] * len(full_dims)) + ind_red2 = np.ravel_multi_index(ind2, [2] * len(full_dims)) + ind_full1 = np.ravel_multi_index(ind1, full_dims) + ind_full2 = np.ravel_multi_index(ind2, full_dims) + + if state.isoper: + assert( + reduced_state[ind_red1, ind_red2] == + state[ind_full1, ind_full2] + ) + elif state.isket: + assert( + reduced_state[ind_red1, 0] == + state[ind_full1, 0] + ) + elif state.isbra: + assert( + reduced_state[0, ind_red2] == + state[0, ind_full2] + ) + + @pytest.mark.parametrize( + "state, full_dims", + [(qutip.rand_dm(8, dims=[[2, 2, 2], [2, 2, 2]]), [3, 2, 3]), + (qutip.rand_ket(8, dims=[[2, 2, 2], [1, 1, 1]]), [2, 3, 3]), + (qutip.Qobj( + qutip.rand_ket(8).full().transpose(), + dims=[[1, 1, 1], [2, 2, 2]]), + [3, 2, 3] + ) + ] + ) + def test_state_expansion(self, state, full_dims): + full_state = expand_qubit_state(state, full_dims) + for _ in range(5): + ind_red1, ind_red2 = \ + np.random.randint(2 ** len(full_dims), size=2) + reduced_dims = [2] * len(full_dims) + ind_full1 = np.ravel_multi_index( + np.unravel_index(ind_red1, reduced_dims), full_dims) + ind_full2 = np.ravel_multi_index( + np.unravel_index(ind_red2, reduced_dims), full_dims) -if __name__ == "__main__": - run_module_suite() + if state.isoper: + assert( + state[ind_red1, ind_red2] == + full_state[ind_full1, ind_full2] + ) + elif state.isket: + assert( + state[ind_red1, 0] == + full_state[ind_full1, 0] + ) + elif state.isbra: + assert( + state[0, ind_red2] == + full_state[0, ind_full2] + ) From e2a89ed880e7db3b39f88a52e21fdbb98c99fb0b Mon Sep 17 00:00:00 2001 From: Boxi Li Date: Sun, 22 Aug 2021 23:31:41 +0200 Subject: [PATCH 2/2] improve the qubit_states function Support a more general way of definition: - use a list or a string for zero/one/plus/minus state - use a list of integers for defining arbitrary disentangled quantum states (up to a global phase). --- src/qutip_qip/qubits.py | 92 ++++++++++++++++++++++++++------- tests/test_optpulseprocessor.py | 8 +-- tests/test_processor.py | 10 ++-- tests/test_qubits.py | 26 +++++----- 4 files changed, 96 insertions(+), 40 deletions(-) diff --git a/src/qutip_qip/qubits.py b/src/qutip_qip/qubits.py index 70568f00c..1e70364cb 100644 --- a/src/qutip_qip/qubits.py +++ b/src/qutip_qip/qubits.py @@ -35,37 +35,93 @@ 'qubit_states', 'truncate_to_qubit_state', 'expand_qubit_state'] import qutip -from qutip import (tensor, basis) import numpy as np from numpy import sqrt -def qubit_states(N=1, states=[0]): +def qubit_states(states): """ - Function to define initial state of the qubits. + Shortcut to generate disentangled qubit states. Parameters ---------- - N : Integer - Number of qubits in the register. - states : List - Initial state of each qubit. + states : list or str + - If a list consisting of ``0``, ``1``, ``"0"``, ``"1"``, ``"+"`` + and ``"-"``, return the corresponding zero/one/plus/minus state. + - If a string consisting of ``0``, ``1``, ``+``, ``-``, + same as above. + - If a list of float or complex numbers, + each number is mapped to a state of the form + :math:`\\sqrt{1 - |a|^2} \\left|0\\right\\rangle + a |1\\rangle`, + where :math:`a` is the given number. Returns - ---------- - qstates : Qobj - List of qubits. + ------- + quantum_states : :obj:`qutip.Qobj` + The generated qubit states. + Examples + -------- + >>> from qutip_qip.qubits import qubit_states + >>> qubit_states([0, 0]) # doctest: +NORMALIZE_WHITESPACE + Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket + Qobj data = + [[1.] + [0.] + [0.] + [0.]] + >>> qubit_states([1, "+"]) # doctest: +NORMALIZE_WHITESPACE + Quantum object: dims = [[2, 2], [1, 1]], shape = + (4, 1), type = ket + Qobj data = + [[0. ] + [0. ] + [0.70710678] + [0.70710678]] + >>> qubit_states("-") # doctest: +NORMALIZE_WHITESPACE + Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket + Qobj data = + [[ 0.70710678] + [-0.70710678]] + >>> qubit_states("1-") # doctest: +NORMALIZE_WHITESPACE + Quantum object: dims = [[2, 2], [1, 1]], shape = + (4, 1), type = ket + Qobj data = + [[ 0. ] + [ 0. ] + [ 0.70710678] + [-0.70710678]] + >>> import numpy as np + >>> qubit_states([1.j/np.sqrt(2)]) # doctest: +NORMALIZE_WHITESPACE + Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket + Qobj data = + [[0.70710678+0.j ] + [0. +0.70710678j]] """ - state_list = [] - for i in range(N): - if N > len(states) and i >= len(states): - state_list.append(0) - else: - state_list.append(states[i]) + if all([np.issubdtype(type(s), np.integer) for s in states]): + return qutip.basis([2] * len(states), states) + + states_map = { + 0: qutip.basis(2, 0), + 1: qutip.basis(2, 1), + "0": qutip.basis(2, 0), + "1": qutip.basis(2, 1), + "+": (qutip.basis(2, 0) + qutip.basis(2, 1)).unit(), + "-": (qutip.basis(2, 0) - qutip.basis(2, 1)).unit(), + } - return tensor([alpha * basis(2, 1) + sqrt(1 - alpha**2) * basis(2, 0) - for alpha in state_list]) + states_list = [] + for s in states: + if s in states_map: + states_list.append(states_map[s]) + elif np.isscalar(s) and abs(s) <= 1: + states_list.append( + s * qutip.basis(2, 1) + + np.sqrt(1 - abs(s)**2) * qutip.basis(2, 0) + ) + else: + raise TypeError(f"Invalid input {s}.") + return qutip.tensor(states_list) def _find_reduced_indices(dims): diff --git a/tests/test_optpulseprocessor.py b/tests/test_optpulseprocessor.py index 4e8105191..c0480fa07 100644 --- a/tests/test_optpulseprocessor.py +++ b/tests/test_optpulseprocessor.py @@ -66,8 +66,8 @@ def test_simple_hadamard(self): qc, num_tslots=num_tslots, evo_time=evo_time, verbose=True) # test run_state - rho0 = qubit_states(1, [0]) - plus = (qubit_states(1, [0]) + qubit_states(1, [1])).unit() + rho0 = qubit_states([0]) + plus = (qubit_states([0]) + qubit_states([1])).unit() result = test.run_state(rho0) assert_allclose(fidelity(result.states[-1], plus), 1, rtol=1.0e-6) @@ -112,8 +112,8 @@ def test_multi_qubits(self): qc = [tensor([identity(2), cnot()])] test.load_circuit(qc, num_tslots=num_tslots, evo_time=evo_time, min_fid_err=1.0e-6) - rho0 = qubit_states(3, [1, 1, 1]) - rho1 = qubit_states(3, [1, 1, 0]) + rho0 = qubit_states([1, 1, 1]) + rho1 = qubit_states([1, 1, 0]) result = test.run_state( rho0, options=Options(store_states=True)) assert_(fidelity(result.states[-1], rho1) > 1-1.0e-6) diff --git a/tests/test_processor.py b/tests/test_processor.py index 4e0c91576..2dd01adac 100644 --- a/tests/test_processor.py +++ b/tests/test_processor.py @@ -278,7 +278,7 @@ def testNoise(self): Test for Processor with noise """ # setup and fidelity without noise - init_state = qubit_states(2, [0, 0, 0, 0]) + init_state = qubit_states([0, 0]) tlist = np.array([0., np.pi/2.]) a = destroy(2) proc = Processor(N=2) @@ -287,7 +287,7 @@ def testNoise(self): proc.pulses[0].coeff = np.array([1]) result = proc.run_state(init_state=init_state) assert_allclose( - fidelity(result.states[-1], qubit_states(2, [0, 1, 0, 0])), + fidelity(result.states[-1], qubit_states([0, 1])), 1, rtol=1.e-7) # decoherence noise @@ -295,7 +295,7 @@ def testNoise(self): proc.add_noise(dec_noise) result = proc.run_state(init_state=init_state) assert_allclose( - fidelity(result.states[-1], qubit_states(2, [0, 1, 0, 0])), + fidelity(result.states[-1], qubit_states([0, 1])), 0.981852, rtol=1.e-3) # white random noise @@ -336,7 +336,7 @@ def testDrift(self): def testChooseSolver(self): # setup and fidelity without noise - init_state = qubit_states(2, [0, 0, 0, 0]) + init_state = qubit_states([0, 0]) tlist = np.array([0., np.pi/2.]) a = destroy(2) proc = Processor(N=2) @@ -345,7 +345,7 @@ def testChooseSolver(self): proc.pulses[0].coeff = np.array([1]) result = proc.run_state(init_state=init_state, solver="mcsolve") assert_allclose( - fidelity(result.states[-1], qubit_states(2, [0, 1, 0, 0])), + fidelity(result.states[-1], qubit_states([0, 1])), 1, rtol=1.e-7) def test_no_saving_intermidiate_state(self): diff --git a/tests/test_qubits.py b/tests/test_qubits.py index 5e6250dff..f7c459659 100644 --- a/tests/test_qubits.py +++ b/tests/test_qubits.py @@ -30,11 +30,10 @@ # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ############################################################################### -from numpy.testing import assert_, run_module_suite import numpy as np import pytest import qutip -from qutip import (tensor, basis) +from qutip import tensor, basis from qutip_qip.qubits import ( qubit_states, truncate_to_qubit_state, expand_qubit_state) @@ -47,17 +46,18 @@ def testQubitStates(self): """ Tests the qubit_states function. """ - psi0_a = basis(2, 0) - psi0_b = qubit_states() - assert_(psi0_a == psi0_b) - - psi1_a = basis(2, 1) - psi1_b = qubit_states(states=[1]) - assert_(psi1_a == psi1_b) - - psi01_a = tensor(psi0_a, psi1_a) - psi01_b = qubit_states(N=2, states=[0, 1]) - assert_(psi01_a == psi01_b) + assert(qubit_states([0]) == basis(2, 0)) + assert(qubit_states([1]) == basis(2, 1)) + assert(qubit_states([0, 1]) == tensor(basis(2, 0), basis(2, 1))) + plus = (basis(2, 0) + basis(2, 1)).unit() + minus = (basis(2, 0) - basis(2, 1)).unit() + assert(qubit_states("-+") == tensor(minus, plus)) + assert(qubit_states("0+") == tensor(basis(2, 0), plus)) + assert(qubit_states("+11") == tensor(plus, basis(2, 1), basis(2, 1))) + assert( + qubit_states([1.j/np.sqrt(2), 1.]) == + tensor(qutip.Qobj([[1/np.sqrt(2)], [1.j/np.sqrt(2)]]), basis(2, 1)) + ) @pytest.mark.parametrize( "state, full_dims",