# This file is part of QuTiP: Quantum Toolbox in Python.
#
# Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# 1. Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
# of its contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
###############################################################################
"""
This module is a collection of random state and operator generators.
The sparsity of the ouput Qobj's is controlled by varing the
`density` parameter.
"""
__all__ = ['rand_herm', 'rand_unitary', 'rand_ket', 'rand_dm']
from scipy import arcsin, sqrt, pi
import numpy as np
import scipy.sparse as sp
from qutip.qobj import Qobj
from qutip.operators import create, destroy, jmat
[docs]def rand_herm(N, density=0.75, dims=None):
"""Creates a random NxN sparse Hermitian quantum object.
Uses :math:`H=X+X^{+}` where :math:`X` is
a randomly generated quantum operator with a given `density`.
Parameters
----------
N : int
Shape of output quantum operator.
density : float
Density between [0,1] of output Hermitian operator.
dims : list
Dimensions of quantum object. Used for specifying
tensor structure. Default is dims=[[N],[N]].
Returns
-------
oper : qobj
NxN Hermitian quantum operator.
"""
if dims:
_check_dims(dims, N, N)
# to get appropriate density of output
# Hermitian operator must convert via:
herm_density = 2.0 * arcsin(density) / pi
X = sp.rand(N, N, herm_density, format='csr')
X.data = X.data - 0.5
Y = X.copy()
Y.data = 1.0j * np.random.random(len(X.data)) - (0.5 + 0.5j)
X = X + Y
X.sort_indices()
X = Qobj(X)
if dims:
return Qobj((X + X.dag()) / 2.0, dims=dims, shape=[N, N])
else:
return Qobj((X + X.dag()) / 2.0)
[docs]def rand_unitary(N, density=0.75, dims=None):
"""Creates a random NxN sparse unitary quantum object.
Uses :math:`\exp(-iH)` where H is a randomly generated
Hermitian operator.
Parameters
----------
N : int
Shape of output quantum operator.
density : float
Density between [0,1] of output Unitary operator.
dims : list
Dimensions of quantum object. Used for specifying
tensor structure. Default is dims=[[N],[N]].
Returns
-------
oper : qobj
NxN Unitary quantum operator.
"""
if dims:
_check_dims(dims, N, N)
U = (-1.0j * rand_herm(N, density)).expm()
U.data.sort_indices()
if dims:
return Qobj(U, dims=dims, shape=[N, N])
else:
return Qobj(U)
[docs]def rand_ket(N, density=1, dims=None):
"""Creates a random Nx1 sparse ket vector.
Parameters
----------
N : int
Number of rows for output quantum operator.
density : float
Density between [0,1] of output ket state.
dims : list
Dimensions of quantum object. Used for specifying
tensor structure. Default is dims=[[N],[1]].
Returns
-------
oper : qobj
Nx1 ket state quantum operator.
"""
if dims:
_check_dims(dims, N, 1)
X = sp.rand(N, 1, density, format='csr')
X.data = X.data - 0.5
Y = X.copy()
Y.data = 1.0j * np.random.random(len(X.data)) - (0.5 + 0.5j)
X = X + Y
X.sort_indices()
X = Qobj(X)
if dims:
return Qobj(X / X.norm(), dims=dims, shape=[N, 1])
else:
return Qobj(X / X.norm())
[docs]def rand_dm(N, density=0.75, pure=False, dims=None):
"""Creates a random NxN density matrix.
Parameters
----------
N : int
Shape of output density matrix.
density : float
Density between [0,1] of output density matrix.
dims : list
Dimensions of quantum object. Used for specifying
tensor structure. Default is dims=[[N],[N]].
Returns
-------
oper : qobj
NxN density matrix quantum operator.
Notes
-----
For small density matrices., choosing a low density will result in an error
as no diagonal elements will be generated such that :math:`Tr(\\rho)=1`.
"""
if dims:
_check_dims(dims, N, N)
if pure:
dm_density = sqrt(density)
psi = rand_ket(N, dm_density)
H = psi * psi.dag()
else:
density = density ** 2
non_zero = 0
tries = 0
while non_zero == 0 and tries < 10:
H = rand_herm(N, density)
H = H.dag() * H
non_zero = sum([H.tr()])
tries += 1
if tries >= 10:
raise ValueError(
"Requested density is too low to generate density matrix.")
H.data.sort_indices()
if dims:
return Qobj(H / H.tr(), dims=dims, shape=[N, N])
else:
return Qobj(H / H.tr())
def rand_kraus_map(N, dims=None):
"""
Creates a random CPTP map on an N-dimensional Hilbert space in Kraus
form.
Parameters
----------
N : int
Length of input/output density matrix.
dims : list
Dimensions of quantum object. Used for specifying
tensor structure. Default is dims=[[N],[N]].
Returns
-------
oper_list : list of qobj
N^2 x N x N qobj operators.
"""
if dims:
_check_dims(dims, N, N)
# Random unitary (Stinespring Dilation)
big_unitary = rand_unitary(N ** 3).data.todense()
orthog_cols = np.array(big_unitary[:, :N])
oper_list = np.reshape(orthog_cols, (N ** 2, N, N))
return list(map(lambda x: Qobj(inpt=x, dims=dims), oper_list))
def rand_super(dim=5):
H = rand_herm(dim)
return propagator(H, np.random.rand(), [
create(dim), destroy(dim), jmat(float(dim - 1) / 2.0, 'z')
])
def _check_dims(dims, N1, N2):
if len(dims) != 2:
raise Exception("Qobj dimensions must be list of length 2.")
if (not isinstance(dims[0], list)) or (not isinstance(dims[1], list)):
raise TypeError(
"Qobj dimension components must be lists. i.e. dims=[[N],[N]]")
if np.prod(dims[0]) != N1 or np.prod(dims[1]) != N2:
raise ValueError("Qobj dimensions must match matrix shape.")
if len(dims[0]) != len(dims[1]):
raise TypeError("Qobj dimension components must have same length.")
# TRAILING IMPORTS
# qutip.propagator depends on rand_dm, so we need to put this import last.
from qutip.propagator import propagator