diff --git a/qutip/core/__init__.py b/qutip/core/__init__.py index 1a7707c6a8..1a33f93297 100644 --- a/qutip/core/__init__.py +++ b/qutip/core/__init__.py @@ -14,6 +14,7 @@ from .blochredfield import * from .energy_restricted import * from .properties import * +from .transformation import * from . import gates del cy # File in cy are not public facing diff --git a/qutip/core/cy/_element.pxd b/qutip/core/cy/_element.pxd index 745df68e9d..4af82c5042 100644 --- a/qutip/core/cy/_element.pxd +++ b/qutip/core/cy/_element.pxd @@ -17,11 +17,23 @@ cdef class _ConstantElement(_BaseElement): cdef readonly object _qobj +cdef class _ConstantSuperElement(_BaseElement): + cdef readonly object _pre + cdef readonly object _post + cdef readonly object _pre_data + cdef readonly object _post_data + cdef readonly object _qobj + + cdef class _EvoElement(_BaseElement): cdef readonly object _qobj cdef readonly Coefficient _coefficient +cdef class _EvoSuperElement(_ConstantSuperElement): + cdef readonly Coefficient _coefficient + + cdef class _FuncElement(_BaseElement): cdef readonly object _func cdef readonly dict _args diff --git a/qutip/core/cy/_element.pyx b/qutip/core/cy/_element.pyx index d77f3e6dfd..e457385032 100644 --- a/qutip/core/cy/_element.pyx +++ b/qutip/core/cy/_element.pyx @@ -5,16 +5,22 @@ #cython: cdvision=True #cython: c_api_binop_methods=True +import warnings from .. import data as _data from qutip.core.cy.coefficient import coefficient_function_parameters +from qutip.core.dimensions import Dimensions +from qutip.core.qobj import Qobj from qutip.core.data cimport Dense, Data, dense from qutip.core.data.matmul cimport * from math import nan as Nan cdef extern from "" namespace "std" nogil: double complex conj(double complex x) -__all__ = ['_ConstantElement', '_EvoElement', - '_FuncElement', '_MapElement', '_ProdElement'] +__all__ = [ + '_ConstantElement', '_EvoElement', + '_ConstantSuperElement', '_EvoSuperElement', + '_FuncElement', '_MapElement', '_ProdElement' +] cdef class _BaseElement: @@ -321,6 +327,130 @@ cdef class _ConstantElement(_BaseElement): def dtype(self): return type(self._data) +cdef class _ConstantSuperElement(_BaseElement): + """ + Constant part of a list format :obj:`.QobjEvo`. + A constant :obj:`.QobjEvo` will contain one `_ConstantElement`:: + + qevo = QobjEvo(H0) + qevo.elements = [_ConstantElement(H0)] + """ + def __init__(self, qobjs): + self._pre = qobjs[0] + self._post = qobjs[1] + self._pre_data = self._pre.data + self._post_data = self._post.data + self._data = None + self._qobj = None + + def __mul__(left, right): + cdef _ConstantSuperElement base + cdef object factor + if type(left) is _ConstantSuperElement: + base = left + factor = right + if type(right) is _ConstantSuperElement: + base = right + factor = left + return _ConstantSuperElement((base._pre * right, base._post)) + + def __matmul__(left, right): + if ( + type(left) is _ConstantSuperElement + and type(right) is _ConstantSuperElement + ): + pre_left = (<_ConstantSuperElement> left)._pre + pre_right = (<_ConstantSuperElement> right)._pre + post_left = (<_ConstantSuperElement> left)._post + post_right = (<_ConstantSuperElement> right)._post + return _ConstantSuperElement( + (pre_left @ pre_right, post_right @ post_left) + ) + elif ( + type(left) is _EvoSuperElement + or type(right) is _EvoSuperElement + ): + return NotImplemented + elif type(left) is _ConstantSuperElement: + return _ConstantElement(left.qobj(0.)) @ right + elif type(right) is _ConstantSuperElement: + return left @ _ConstantElement(right.qobj(0.)) + raise Exception("Should never reach this.") + + cpdef Data data(self, t): + if self._data is None: + self._data = _data.kron_transpose(self._post_data, self._pre_data) + return self._data + + cpdef object qobj(self, t): + if self._qobj is None: + dims = Dimensions.from_prepost(self._pre._dims, self._post._dims) + self._qobj = Qobj( + self.data(t), + dims=dims, + superrep='super', + isherm=self._pre._isherm and self._post._isherm, + copy=False + ) + return self._qobj + + cpdef object coeff(self, t): + return 1. + + cdef Data matmul_data_t(_ConstantSuperElement self, t, Data state, Data out=None): + if state.shape[1] == 1 and self._post_data.shape[0] != 1: + # Operator ket: restack + if type(state) == Dense: + state = _data.column_unstack_dense( + state, self._pre_data.shape[1], inplace=state.fortran + ) + else: + state = _data.column_unstack(state, self._pre_data.shape[1]) + if type(out) is Dense: + out = _data.column_unstack_dense( + out, self._pre_data.shape[1], inplace=out.fortran + ) + elif out is not None: + out = _data.column_unstack(out, self._pre_data.shape[1]) + + try: + out = self.matmul_data_t(t, state, out) + finally: + if state.fortran: + # `state` was reshaped inplace, restore it's original shape + _data.column_stack_dense(state, inplace=state.fortran) + if type(out) is Dense: + out = _data.column_stack_dense(out, inplace=out.fortran) + else: + out = _data.column_stack(out) + return out + + # Core computation of matmul + rho_post = _data.matmul(state, self._post_data) + if out is None: + return _data.matmul(self._pre_data, rho_post) + elif type(state) is Dense and type(out) is Dense: + imatmul_data_dense(self._pre_data, rho_post, 1., out=out) + return out + else: + return _data.add(out, _data.matmul(self._pre_data, rho_post)) + + def linear_map(self, f, anti=False): + warnings.warn("Conversion from spre / spost to super", UserWarning) + return _ConstantElement(f(self.qobj(0.))) + + def replace_arguments(self, args, cache=None): + return self + + def __call__(self, t, args=None): + return self.qobj(0.) + + @property + def dtype(self): + if type(self._pre_data) is type(self._post_data): + return type(self._pre) + return None + cdef class _EvoElement(_BaseElement): """ @@ -383,6 +513,120 @@ cdef class _EvoElement(_BaseElement): return type(self._data) +cdef class _EvoSuperElement(_ConstantSuperElement): + """ + A pair of a :obj:`.Qobj` and a :obj:`.Coefficient` from the list format + time-dependent operator:: + + qevo = QobjEvo([[H0, coeff0], [H1, coeff1]]) + qevo.elements = [_EvoElement(H0, coeff0), _EvoElement(H1, coeff1)] + """ + def __init__(self, qobjs, coefficient): + self._pre = qobjs[0] + self._post = qobjs[1] + self._pre_data = self._pre.data + self._post_data = self._post.data + self._data = None + self._qobj = None + self._coefficient = coefficient + + def __mul__(left, right): + cdef _EvoSuperElement base + cdef object factor + if type(left) is _EvoSuperElement: + base = left + factor = right + if type(right) is _EvoSuperElement: + base = right + factor = left + return _EvoSuperElement( + (base._pre * factor, base._post), + base._coefficient + ) + + def __matmul__(left, right): + if not isinstance(right, _ConstantSuperElement): + return left @ _EvoElement(right.qobj(0.), right._ceofficient) + if not isinstance(left, _ConstantSuperElement): + return _EvoElement(left.qobj(0.), left._ceofficient) @ right + + pre_left = (<_ConstantSuperElement> left)._pre + pre_right = (<_ConstantSuperElement> right)._pre + post_left = (<_ConstantSuperElement> left)._post + post_right = (<_ConstantSuperElement> right)._post + + if isinstance(left, _EvoSuperElement) and isinstance(right, _EvoSuperElement): + coefficient = left._coefficient * right._coefficient + elif isinstance(left, _EvoSuperElement): + coefficient = left._coefficient + elif isinstance(right, _EvoSuperElement): + coefficient = right._coefficient + else: + raise Exception("Should not be here!") + + return _EvoSuperElement( + (pre_left @ pre_right, post_right @ post_left), + coefficient + ) + + cpdef object coeff(self, t): + return self._coefficient(t) + + cdef Data matmul_data_t(_EvoSuperElement self, t, Data state, Data out=None): + if state.shape[1] == 1 and self._post_data.shape[0] != 1: + # Operator ket: restack + if type(state) == Dense: + # 1D matrix are both C and fortran. + # So we can always stack them inplace. + ( state).fortran = True + state = _data.column_unstack_dense( + state, self._pre_data.shape[1], inplace=True + ) + else: + state = _data.column_unstack(state, self._pre_data.shape[1]) + if type(out) is Dense: + ( out).fortran = True + out = _data.column_unstack_dense( + out, self._pre_data.shape[1], inplace=True + ) + elif out is not None: + out = _data.column_unstack(out, self._pre_data.shape[1]) + + try: + out = self.matmul_data_t(t, state, out) + finally: + # `state` was reshaped inplace, restore it's original shape + _data.column_stack_dense(state, inplace=True) + if type(out) is Dense: + out = _data.column_stack_dense(out, inplace=True) + else: + out = _data.column_stack(out) + return out + + # Core computation of matmul + rho_post = _data.matmul(state, self._post_data, self.coeff(t)) + if out is None: + return _data.matmul(self._pre_data, rho_post) + elif type(rho_post) is Dense and type(out) is Dense: + imatmul_data_dense(self._pre_data, rho_post, 1., out=out) + return out + else: + return _data.add(out, _data.matmul(self._pre_data, rho_post)) + + def linear_map(self, f, anti=False): + warnings.warn("Conversion from spre / spost to super", UserWarning) + return _EvoElement( + f(self.qobj(0.)), + self._coefficient.conj() if anti else self._coefficient, + ) + + def replace_arguments(self, args, cache=None): + return _EvoSuperElement( + (self._pre, self._post), + self._coefficient.replace_arguments(args) + ) + + cdef class _FuncElement(_BaseElement): """ Used with :obj:`.QobjEvo` to build an evolution term from a function with diff --git a/qutip/core/cy/qobjevo.pyx b/qutip/core/cy/qobjevo.pyx index a65a99ebad..5c973222eb 100644 --- a/qutip/core/cy/qobjevo.pyx +++ b/qutip/core/cy/qobjevo.pyx @@ -265,20 +265,33 @@ cdef class QobjEvo: """ Read a Q_object item and return an element for that item. """ if isinstance(op, Qobj): out = _ConstantElement(op.copy() if copy else op) - qobj = op - elif isinstance(op, list): + dims = op._dims + elif isinstance(op, tuple) and len(op) == 2: + out = _ConstantSuperElement( + (op[0].copy(), op[1].copy()) if copy else op + ) + dims = Dimensions.from_prepost(op[0]._dims, op[1]._dims) + elif isinstance(op, list) and isinstance(op[0], Qobj): out = _EvoElement( op[0].copy() if copy else op[0], coefficient(op[1], tlist=tlist, args=args, order=order, boundary_conditions=boundary_conditions) ) - qobj = op[0] + dims = op[0]._dims + elif isinstance(op, list) and isinstance(op[0], tuple): + out = _EvoSuperElement( + (op[0][0].copy(), op[0][1].copy()) if copy else op[0], + coefficient(op[1], tlist=tlist, args=args, order=order, + boundary_conditions=boundary_conditions) + ) + dims = Dimensions.from_prepost(op[0][0]._dims, op[0][1]._dims) elif isinstance(op, _BaseElement): out = op - qobj = op.qobj(0) + dims = op.qobj(0)._dims elif callable(op): out = _FuncElement(op, args, style=function_style) qobj = out.qobj(0) + dims = qobj._dims if not isinstance(qobj, Qobj): raise TypeError( "Function based time-dependent elements must have the" @@ -293,12 +306,12 @@ cdef class QobjEvo: ) if self._dims is None: - self._dims = qobj._dims - self.shape = qobj.shape - elif self._dims != qobj._dims: + self._dims = dims + self.shape = dims.shape + elif self._dims != dims: raise ValueError( - f"QobjEvo term {op!r} has dims {qobj.dims!r} and shape" - f" {qobj.shape!r} but previous terms had dims {self.dims!r}" + f"QobjEvo term {op!r} has dims {dims!r} and shape" + f" {dims.shape!r} but previous terms had dims {self.dims!r}" f" and shape {self.shape!r}." ) @@ -547,11 +560,16 @@ cdef class QobjEvo: raise TypeError("incompatible dimensions" + str(self.dims) + ", " + str(other.dims)) self.elements.append(_ConstantElement(other)) + + elif other == 0.: + pass + elif ( isinstance(other, numbers.Number) and self._dims[0] == self._dims[1] ): self.elements.append(_ConstantElement(other * qutip.qeye_like(self))) + else: return NotImplemented return self diff --git a/qutip/core/data/matmul.pyx b/qutip/core/data/matmul.pyx index efaa1047a6..82db10a4cb 100644 --- a/qutip/core/data/matmul.pyx +++ b/qutip/core/data/matmul.pyx @@ -33,6 +33,10 @@ cnp.import_array() cdef extern from *: void *PyMem_Calloc(size_t n, size_t elsize) + +cdef extern from "" namespace "std" nogil: + double complex conj "conj"(double complex x) + # This function is templated over integral types on import to allow `idxint` to # be any signed integer (though likely things will only work for >=32-bit). To # change integral types, you only need to change the `idxint` definitions in @@ -56,6 +60,7 @@ __all__ = [ 'matmul', 'matmul_csr', 'matmul_dense', 'matmul_dia', 'matmul_csr_dense_dense', 'matmul_dia_dense_dense', 'matmul_dense_dia_dense', 'multiply', 'multiply_csr', 'multiply_dense', 'multiply_dia', + 'matmul_dag', 'matmul_dag_data', 'matmul_dag_dense_csr_dense', ] @@ -518,6 +523,79 @@ cpdef Dense matmul_dense_dia_dense(Dense left, Dia right, double complex scale=1 return out +cpdef Dense matmul_dag_dense_csr_dense( + Dense left, CSR right, + double complex scale=1, Dense out=None +): + """ + Perform the operation + ``out := scale * (left @ right) + out`` + where `left`, `right` and `out` are matrices. `scale` is a complex scalar, + defaulting to 1. + + If `out` is not given, it will be allocated as if it were a zero matrix. + """ + if left.shape[1] != right.shape[1]: + raise ValueError( + "incompatible matrix shapes " + + str(left.shape) + + " and " + + str(right.shape) + ) + if ( + out is not None and ( + out.shape[0] != left.shape[0] + or out.shape[1] != right.shape[0] + ) + ): + raise ValueError( + "incompatible output shape, got " + + str(out.shape) + + " but needed " + + str((left.shape[0], right.shape[0])) + ) + cdef Dense tmp = None + if out is None: + out = dense.zeros(left.shape[0], right.shape[0], left.fortran) + if bool(left.fortran) != bool(out.fortran): + msg = ( + "out matrix is {}-ordered".format('Fortran' if out.fortran else 'C') + + " but input is {}-ordered".format('Fortran' if left.fortran else 'C') + ) + warnings.warn(msg, OrderEfficiencyWarning) + # Rather than making loads of copies of the same code, we just moan at + # the user and then transpose one of the arrays. We prefer to have + # `right` in Fortran-order for cache efficiency. + if left.fortran: + tmp = out + out = out.reorder() + else: + left = left.reorder() + cdef idxint row, col, ptr, idx_l, idx_out, out_row + cdef idxint stride_in_col, stride_in_row, stride_out_row, stride_out_col + cdef idxint nrows=left.shape[0], ncols=right.shape[1] + cdef double complex val + stride_in_col = left.shape[0] if left.fortran else 1 + stride_in_row = 1 if left.fortran else left.shape[1] + stride_out_col = out.shape[0] if out.fortran else 1 + stride_out_row = 1 if out.fortran else out.shape[1] + + # A @ B.dag = (B* @ A.T).T + # Todo: make a conj version of _matmul_csr_vector? + for row in range(right.shape[0]): + for ptr in range(right.row_index[row], right.row_index[row + 1]): + val = scale * conj(right.data[ptr]) + col = right.col_index[ptr] + for out_row in range(out.shape[0]): + idx_out = row * stride_out_col + out_row * stride_out_row + idx_l = col * stride_in_col + out_row * stride_in_row + out.data[idx_out] += val * left.data[idx_l] + if tmp is None: + return out + memcpy(tmp.data, out.data, ncols * nrows * sizeof(double complex)) + return tmp + + cpdef CSR multiply_csr(CSR left, CSR right): """Element-wise multiplication of CSR matrices.""" if left.shape[0] != right.shape[0] or left.shape[1] != right.shape[1]: @@ -732,6 +810,48 @@ multiply.add_specialisations([ ], _defer=True) +cpdef Data matmul_dag_data( + Data left, Data right, + double complex scale=1, Dense out=None +): + return matmul(left, right.adjoint(), scale, out) + + +matmul_dag = _Dispatcher( + _inspect.Signature([ + _inspect.Parameter('left', _inspect.Parameter.POSITIONAL_ONLY), + _inspect.Parameter('right', _inspect.Parameter.POSITIONAL_ONLY), + _inspect.Parameter('scale', _inspect.Parameter.POSITIONAL_OR_KEYWORD, + default=1), + ]), + name='matmul_dag', + module=__name__, + inputs=('left', 'right'), + out=True, +) +matmul_dag.__doc__ =\ + """ + Compute the matrix multiplication of two matrices, with the operation + scale * (left @ right.dag) + where `scale` is (optionally) a scalar, and `left` and `right` are + matrices. + + Parameters + ---------- + left : Data + The left operand as either a bra or a ket matrix. + + right : Data + The right operand as a ket matrix. + + scale : complex, optional + The scalar to multiply the output by. + """ +matmul_dag.add_specialisations([ + (Dense, CSR, Dense, matmul_dag_dense_csr_dense), + (Data, Data, Data, matmul_dag_data), +], _defer=True) + del _inspect, _Dispatcher diff --git a/qutip/core/data/properties.pyx b/qutip/core/data/properties.pyx index 2ab25ca543..37344925be 100644 --- a/qutip/core/data/properties.pyx +++ b/qutip/core/data/properties.pyx @@ -339,6 +339,11 @@ cpdef bint isequal_dia(Dia A, Dia B, double atol=-1, double rtol=-1): cdef double complex *ptr_b cdef idxint size=A.shape[1] + if A.num_diag == 0: + return iszero_dia(B, atol) + if B.num_diag == 0: + return iszero_dia(A, atol) + # TODO: # Works only for a sorted offsets list. # We don't have a check for whether it's already sorted, but it should be diff --git a/qutip/core/dimensions.py b/qutip/core/dimensions.py index 56fe353964..3065199c20 100644 --- a/qutip/core/dimensions.py +++ b/qutip/core/dimensions.py @@ -835,6 +835,18 @@ def __init__(self, from_: Space, to_: Space): self.superrep = 'mixed' self.__setitem__ = _frozen + @staticmethod + def from_prepost(pre, post): + pre = Dimensions(pre) + post = Dimensions(post) + return Dimensions( + SuperSpace(Dimensions(pre[1], post[0])), + SuperSpace(Dimensions(pre[0], post[1])), + ) + + def transpose(self): + return Dimensions(self.to_, self.from_) + def __eq__(self, other: "Dimensions") -> bool: if isinstance(other, Dimensions): return ( diff --git a/qutip/core/energy_restricted.py b/qutip/core/energy_restricted.py index 1a12f93b3e..723961f502 100644 --- a/qutip/core/energy_restricted.py +++ b/qutip/core/energy_restricted.py @@ -182,7 +182,7 @@ def enr_thermal_dm(dims, excitations, n, *, dtype=None): diags /= np.sum(diags) out = qdiags(diags, 0, dims=enr_dims, shape=(nstates, nstates), dtype=dtype) - out._isherm = True + out.isherm = True return out diff --git a/qutip/core/expect.py b/qutip/core/expect.py index 647aba3c69..59623c7a10 100644 --- a/qutip/core/expect.py +++ b/qutip/core/expect.py @@ -78,7 +78,7 @@ def _single_qobj_expect(oper, state): Private function used by expect to calculate expectation values of Qobjs. """ if not oper.isoper or not (state.isket or state.isoper): - raise TypeError('invalid operand types') + raise TypeError(f'invalid operand types, {oper.type} and {state.type}') if oper.dims[1] != state.dims[0]: msg = ( "incompatible dimensions " diff --git a/qutip/core/gates.py b/qutip/core/gates.py index 36d15e1074..80b5d0f332 100644 --- a/qutip/core/gates.py +++ b/qutip/core/gates.py @@ -840,5 +840,5 @@ def qubit_clifford_group(*, dtype: LayerType = None) -> list[Qobj]: ] for gate in gates: gate.isherm - gate._isunitary = True + gate.isunitary = True return gates diff --git a/qutip/core/metrics.py b/qutip/core/metrics.py index 011536108b..56dceab2ae 100644 --- a/qutip/core/metrics.py +++ b/qutip/core/metrics.py @@ -109,6 +109,8 @@ def _hilbert_space_dims(oper): return oper.dims elif oper.type == 'super' and oper.superrep in ['choi', 'chi', 'super']: return [oper.dims[0][1], oper.dims[1][0]] + elif oper.type == 'super' and oper.superrep in ['kraus']: + return oper._pre_dims().as_list() else: raise TypeError('oper is not a valid quantum channel!') diff --git a/qutip/core/operators.py b/qutip/core/operators.py index c59765571e..9677ac2325 100644 --- a/qutip/core/operators.py +++ b/qutip/core/operators.py @@ -341,11 +341,11 @@ def spin_J_set(j: float, *, dtype: LayerType = None) -> tuple[Qobj]: _SIGMAP = jmat(0.5, '+') _SIGMAM = jmat(0.5, '-') _SIGMAX = 2 * jmat(0.5, 'x') -_SIGMAX._isunitary = True +_SIGMAX.isunitary = True _SIGMAY = 2 * jmat(0.5, 'y') -_SIGMAY._isunitary = True +_SIGMAY.isunitary = True _SIGMAZ = 2 * jmat(0.5, 'z') -_SIGMAZ._isunitary = True +_SIGMAZ.isunitary = True def sigmap(*, dtype: LayerType = None) -> Qobj: @@ -695,7 +695,7 @@ def _f_op(n_sites, site, action, dtype: LayerType = None,): opers = [s_z] * site + [operator] + [eye] * (n_sites - site - 1) out = tensor(opers).to(dtype) out.isherm = False - out._isunitary = False + out.isunitary = False return out @@ -869,7 +869,7 @@ def position(N: int, offset: int = 0, *, dtype: LayerType = None) -> Qobj: a = destroy(N, offset=offset, dtype=dtype) position = np.sqrt(0.5) * (a + a.dag()) position.isherm = True - position._isunitary = False + position.isunitary = False return position.to(dtype) @@ -899,7 +899,7 @@ def momentum(N: int, offset: int = 0, *, dtype: LayerType = None) -> Qobj: a = destroy(N, offset=offset, dtype=dtype) momentum = -1j * np.sqrt(0.5) * (a - a.dag()) momentum.isherm = True - momentum._isunitary = False + momentum.isunitary = False return momentum.to(dtype) @@ -988,7 +988,7 @@ def squeeze( op = 0.5*np.conj(z)*asq - 0.5*z*asq.dag() out = op.expm(dtype=dtype) out.isherm = (N == 2) or (z == 0.) - out._isunitary = True + out.isunitary = True return out @@ -1067,7 +1067,7 @@ def displace( a = destroy(N, offset=offset) out = (alpha * a.dag() - np.conj(alpha) * a).expm(dtype=dtype) out.isherm = (alpha == 0.) - out._isunitary = True + out.isunitary = True return out @@ -1122,12 +1122,12 @@ def qutrit_ops(*, dtype: LayerType = None) -> list[Qobj]: for i in range(3): op = basis[i] @ basis[i].dag() op.isherm = True - op._isunitary = False + op.isunitary = False out.append(op) for i in range(3): op = basis[i] @ basis[(i+1)%3].dag() op.isherm = False - op._isunitary = False + op.isunitary = False out.append(op) return out @@ -1214,7 +1214,7 @@ def charge( Nmin = -Nmax diag = frac * np.arange(Nmin, Nmax+1, dtype=float) out = qdiags(diag, 0, dtype=dtype) - out._isunitary = (len(diag) <= 2) and np.all(np.abs(diag) == 1.) + out.isunitary = (len(diag) <= 2) and np.all(np.abs(diag) == 1.) return out @@ -1243,7 +1243,7 @@ def tunneling(N: int, m: int = 1, *, dtype: LayerType = None) -> Qobj: diags = [np.ones(N-m, dtype=int), np.ones(N-m, dtype=int)] T = qdiags(diags, [m, -m], dtype=dtype) T.isherm = True - T._isunitary = (m * 2 == N) + T.isunitary = (m * 2 == N) return T diff --git a/qutip/core/qobj/__init__.py b/qutip/core/qobj/__init__.py new file mode 100644 index 0000000000..9eb07dcdc6 --- /dev/null +++ b/qutip/core/qobj/__init__.py @@ -0,0 +1,4 @@ +from ._base import Qobj +from ._state import * +from ._operator import * +from ._super import * diff --git a/qutip/core/qobj/_base.py b/qutip/core/qobj/_base.py new file mode 100644 index 0000000000..8ac7242bb3 --- /dev/null +++ b/qutip/core/qobj/_base.py @@ -0,0 +1,936 @@ +from __future__ import annotations + +import numpy as np +import numbers +import functools +import scipy.sparse + +import qutip +from ... import __version__ +from ...settings import settings +from .. import data as _data +from ..dimensions import ( + enumerate_flat, collapse_dims_super, flatten, unflatten, Dimensions +) +from typing import Any, Literal +from qutip.typing import LayerType, DimensionLike +from numpy.typing import ArrayLike + + +def _latex_real(x): + if not x: + return "0" + if not 0.001 <= abs(x) < 1000: + base, exp = "{:.3e}".format(x).split('e') + return base + r"\times10^{{ {:d} }}".format(int(exp)) + if abs(x - int(x)) < 0.001: + return "{:d}".format(round(x)) + return "{:.3f}".format(x) + + +def _latex_complex(x): + if abs(x.imag) < 0.001: + return _latex_real(x.real) + if abs(x.real) < 0.001: + return _latex_real(x.imag) + "j" + sign = "+" if x.imag > 0 else "-" + return "(" + _latex_real(x.real) + sign + _latex_real(abs(x.imag)) + "j)" + + +def _latex_row(row, cols, data): + if row is None: + bits = (r"\ddots" if col is None else r"\vdots" for col in cols) + else: + bits = (r"\cdots" if col is None else _latex_complex(data[row, col]) + for col in cols) + return " & ".join(bits) + + +class _QobjBuilder(type): + qobjtype_to_class = {} + + @staticmethod + def _initialize_data(arg, raw_dims, copy): + flags = {} + if isinstance(arg, _data.Data): + data = arg.copy() if copy else arg + dims = Dimensions(raw_dims or [[arg.shape[0]], [arg.shape[1]]]) + elif isinstance(arg, Qobj): + data = arg.data.copy() if copy else arg.data + dims = Dimensions(raw_dims or arg._dims) + flags["isherm"] = arg._isherm + flags["isunitary"] = arg._isunitary + else: + data = _data.create(arg, copy=copy) + dims = Dimensions( + raw_dims or [[data.shape[0]], [data.shape[1]]] + ) + if dims.shape != data.shape: + raise ValueError('Provided dimensions do not match the data: ' + + f"{dims.shape} vs {data.shape}") + + return data, dims, flags + + def __call__( + cls, + arg: ArrayLike | Any = None, + dims: DimensionLike = None, + copy: bool = True, + superrep: str = None, + isherm: bool = None, + isunitary: bool = None + ): + if cls is not Qobj: + out = cls.__new__(cls) + out.__init__( + arg, dims, + copy=copy, superrep=superrep, + isherm=isherm, isunitary=isunitary + ) + return out + + data, dims, flags = _QobjBuilder._initialize_data(arg, dims, copy) + if isherm is not None: + flags["isherm"] = isherm + if isunitary is not None: + flags["isunitary"] = isunitary + if superrep is not None: + dims = dims.replace_superrep(superrep) + + instance_class = _QobjBuilder.qobjtype_to_class[dims.type] + + new_qobj = instance_class.__new__(instance_class) + new_qobj.__init__(data, dims, **flags) + return new_qobj + + +def _require_equal_type(method): + """ + Decorate a binary Qobj method to ensure both operands are Qobj and of the + same type and dimensions. Promote numeric scalar to identity matrices of + the same type and shape. + """ + @functools.wraps(method) + def out(self, other): + if isinstance(other, Qobj): + if self._dims != other._dims: + msg = ( + "incompatible dimensions " + + repr(self.dims) + " and " + repr(other.dims) + ) + raise ValueError(msg) + return method(self, other) + if other == 0: + return method(self, other) + if self._dims.issquare and isinstance(other, numbers.Number): + scale = complex(other) + other = Qobj(_data.identity(self.shape[0], scale, + dtype=type(self.data)), + dims=self._dims, + isherm=(scale.imag == 0), + isunitary=(abs(abs(scale)-1) < settings.core['atol']), + copy=False) + return method(self, other) + return NotImplemented + + return out + + +class Qobj(metaclass=_QobjBuilder): + _dims: Dimensions + dims: Dimensions + _data: Data + data: data + __array_ufunc__ = None + + + def __init__( + self, + arg: ArrayLike | Any = None, + dims: DimensionLike = None, + copy: bool = True, + superrep: str = None, + isherm: bool = None, + isunitary: bool = None + ): + if not (isinstance(arg, _data.Data) and isinstance(dims, Dimensions)): + arg, dims, flags = _QobjBuilder._initialize_data(arg, dims, copy) + if isherm is None and flags["isherm"] is not None: + isherm = flags["isherm"] + if isunitary is None and flags["isunitary"] is not None: + isunitary = flags["isunitary"] + if superrep is not None: + dims = dims.replace_superrep(superrep) + + self._data = arg + self._dims = dims + self._flags = {} + if isherm is not None: + self._flags["isherm"] = isherm + if isunitary is not None: + self._flags["isunitary"] = isunitary + + @property + def type(self) -> str: + return self._dims.type + + @property + def data(self) -> _data.Data: + return self._data + + @data.setter + def data(self, data: _data.Data): + if not isinstance(data, _data.Data): + raise TypeError('Qobj data must be a data-layer format.') + if self._dims.shape != data.shape: + raise ValueError('Provided data do not match the dimensions: ' + + f"{self._dims.shape} vs {data.shape}") + self._data = data + + @property + def dtype(self): + return type(self._data) + + @property + def dims(self) -> list[list[int]] | list[list[list[int]]]: + return self._dims.as_list() + + @dims.setter + def dims(self, dims: list[list[int]] | list[list[list[int]]] | Dimensions): + dims = Dimensions(dims, rep=self.superrep) + if dims.shape != self._data.shape: + raise ValueError('Provided dimensions do not match the data: ' + + f"{dims.shape} vs {self._data.shape}") + self._dims = dims + + def copy(self) -> Qobj: + """Create identical copy""" + return Qobj(arg=self._data, + dims=self._dims, + isherm=self._isherm, + isunitary=self._isunitary, + copy=True) + + def to(self, data_type: LayerType, copy: bool=False) -> Qobj: + """ + Convert the underlying data store of this `Qobj` into a different + storage representation. + + The different storage representations available are the "data-layer + types" which are known to :obj:`qutip.core.data.to`. By default, these + are :class:`~qutip.core.data.CSR`, :class:`~qutip.core.data.Dense` and + :class:`~qutip.core.data.Dia`, which respectively construct a + compressed sparse row matrix, diagonal matrix and a dense one. Certain + algorithms and operations may be faster or more accurate when using a + more appropriate data store. + + Parameters + ---------- + data_type : type, str + The data-layer type or its string alias that the data of this + :class:`Qobj` should be converted to. + + copy : Bool + If the data store is already in the format requested, whether the + function should return returns `self` or a copy. + + Returns + ------- + Qobj + A :class:`Qobj` with the data stored in the requested format. + """ + data_type = _data.to.parse(data_type) + if type(self._data) is data_type and copy: + return self.copy() + elif type(self._data) is data_type: + return self + return self.__class__( + _data.to(data_type, self._data), + dims=self._dims, + isherm=self._isherm, + isunitary=self._isunitary, + copy=False + ) + + @_require_equal_type + def __add__(self, other: Qobj | complex) -> Qobj: + if other == 0: + return self.copy() + return Qobj(_data.add(self._data, other._data), + dims=self._dims, + isherm=(self._isherm and other._isherm) or None, + copy=False) + + def __radd__(self, other: Qobj | complex) -> Qobj: + return self.__add__(other) + + @_require_equal_type + def __sub__(self, other: Qobj | complex) -> Qobj: + if other == 0: + return self.copy() + return Qobj(_data.sub(self._data, other._data), + dims=self._dims, + isherm=(self._isherm and other._isherm) or None, + copy=False) + + def __rsub__(self, other: Qobj | complex) -> Qobj: + return self.__neg__().__add__(other) + + def __neg__(self) -> Qobj: + return Qobj(_data.neg(self._data), + dims=self._dims, + isherm=self._isherm, + isunitary=self._isunitary, + copy=False) + + def __mul__(self, other: complex) -> Qobj: + """ + If other is a Qobj, we dispatch to __matmul__. If not, we + check that other is a valid complex scalar, i.e., we can do + complex(other). Otherwise, we return NotImplemented. + """ + + if isinstance(other, Qobj): + return self.__matmul__(other) + + # We send other to mul instead of complex(other) to be more flexible. + # The dispatcher can then decide how to handle other and return + # TypeError if it does not know what to do with the type of other. + try: + out = _data.mul(self._data, other) + except TypeError: + return NotImplemented + + # Infer isherm and isunitary if possible + try: + multiplier = complex(other) + isherm = (self._isherm and multiplier.imag == 0) or None + isunitary = (abs(abs(multiplier) - 1) < settings.core['atol'] + if self._isunitary else None) + except TypeError: + isherm = None + isunitary = None + + return Qobj(out, + dims=self._dims, + isherm=isherm, + isunitary=isunitary, + copy=False) + + def __rmul__(self, other: complex) -> Qobj: + # Shouldn't be here unless `other.__mul__` has already been tried, so + # we _shouldn't_ check that `other` is `Qobj`. + return self.__mul__(other) + + def __truediv__(self, other: complex) -> Qobj: + return self.__mul__(1 / other) + + def __matmul__(self, other: Qobj) -> Qobj: + if not isinstance(other, Qobj): + try: + other = Qobj(other) + except TypeError: + return NotImplemented + new_dims = self._dims @ other._dims + if new_dims.type == 'scalar': + return _data.inner(self._data, other._data) + + return Qobj( + _data.matmul(self._data, other._data), + dims=new_dims, + isunitary=self._isunitary and other._isunitary, + copy=False + ) + + def __eq__(self, other) -> bool: + if self is other: + return True + if not isinstance(other, Qobj) or self._dims != other._dims: + return False + # isequal uses both atol and rtol from settings.core + return _data.isequal(self._data, other._data) + + def __and__(self, other: Qobj) -> Qobj: + """ + Syntax shortcut for tensor: + A & B ==> tensor(A, B) + """ + return qutip.tensor(self, other) + + def dag(self) -> Qobj: + """Get the Hermitian adjoint of the quantum object.""" + if self._isherm: + return self.copy() + return Qobj(_data.adjoint(self._data), + dims=Dimensions(self._dims[0], self._dims[1]), + isherm=self._isherm, + isunitary=self._isunitary, + copy=False) + + def conj(self) -> Qobj: + """Get the element-wise conjugation of the quantum object.""" + return Qobj(_data.conj(self._data), + dims=self._dims, + isherm=self._isherm, + isunitary=self._isunitary, + copy=False) + + def trans(self) -> Qobj: + """Get the matrix transpose of the quantum operator. + + Returns + ------- + oper : :class:`.Qobj` + Transpose of input operator. + """ + return Qobj(_data.transpose(self._data), + dims=Dimensions(self._dims[0], self._dims[1]), + isherm=self._isherm, + isunitary=self._isunitary, + copy=False) + + def _str_header(self): + out = ", ".join([ + f"{self.__class__.__name__}: dims={self.dims}", + f"shape={self._data.shape}", + f"type={repr(self.type)}", + f"dtype={self.dtype.__name__}", + ]) + # TODO: Should this be here? + if self.type in ('oper', 'super'): + out += ", isherm=" + str(self.isherm) + if self.issuper and self.superrep != 'super': + out += ", superrep=" + repr(self.superrep) + return out + + def __str__(self): + if self.data.shape[0] * self.data.shape[0] > 100_000_000: + # If the system is huge, don't attempt to convert to a dense matrix + # and then to string, because it is pointless and is likely going + # to produce memory errors. Instead print the sparse data string + # representation. + data = _data.to(_data.CSR, self.data).as_scipy() + elif _data.iszero(_data.sub(self.data.conj(), self.data)): + # TODO: that check could be slow... + data = np.real(self.full()) + else: + data = self.full() + return "\n".join([self._str_header(), "Qobj data =", str(data)]) + + def __repr__(self): + # give complete information on Qobj without print statement in + # command-line we cant realistically serialize a Qobj into a string, + # so we simply return the informal __str__ representation instead.) + return self.__str__() + + def _repr_latex_(self): + """ + Generate a LaTeX representation of the Qobj instance. Can be used for + formatted output in ipython notebook. + """ + half_length = 5 + n_rows, n_cols = self.data.shape + # Choose which rows and columns we're going to output, or None if that + # element should be truncated. + rows = list(range(min((half_length, n_rows)))) + if n_rows <= half_length * 2: + rows += list(range(half_length, min((2*half_length, n_rows)))) + else: + rows.append(None) + rows += list(range(n_rows - half_length, n_rows)) + cols = list(range(min((half_length, n_cols)))) + if n_cols <= half_length * 2: + cols += list(range(half_length, min((2*half_length, n_cols)))) + else: + cols.append(None) + cols += list(range(n_cols - half_length, n_cols)) + # Make the data array. + data = r'$$\left(\begin{array}{cc}' + data += r"\\".join(_latex_row(row, cols, self.data.to_array()) + for row in rows) + data += r'\end{array}\right)$$' + return self._str_header() + data + + def __getstate__(self): + # defines what happens when Qobj object gets pickled + self.__dict__.update({'qutip_version': __version__[:5]}) + return self.__dict__ + + def __setstate__(self, state): + # defines what happens when loading a pickled Qobj + # TODO: what happen with miss matched child class? + if 'qutip_version' in state.keys(): + del state['qutip_version'] + (self.__dict__).update(state) + + def __getitem__(self, ind): + # TODO: should we require that data-layer types implement this? This + # isn't the right way of handling it, for sure. + if isinstance(self._data, _data.CSR): + data = self._data.as_scipy() + elif isinstance(self._data, _data.Dense): + data = self._data.as_ndarray() + else: + data = self._data + try: + out = data[ind] + return out.toarray() if scipy.sparse.issparse(out) else out + except TypeError: + pass + return data.to_array()[ind] + + def full( + self, + order: Literal['C', 'F'] = 'C', + squeeze: bool = False + ) -> np.ndarray: + """Dense array from quantum object. + + Parameters + ---------- + order : str {'C', 'F'} + Return array in C (default) or Fortran ordering. + squeeze : bool {False, True} + Squeeze output array. + + Returns + ------- + data : array + Array of complex data from quantum objects `data` attribute. + """ + out = np.asarray(self.data.to_array(), order=order) + return out.squeeze() if squeeze else out + + def data_as(self, format: str = None, copy: bool = True) -> Any: + """Matrix from quantum object. + + Parameters + ---------- + format : str, default: None + Type of the output, "ndarray" for ``Dense``, "csr_matrix" for + ``CSR``. A ValueError will be raised if the format is not + supported. + + copy : bool {False, True} + Whether to return a copy + + Returns + ------- + data : numpy.ndarray, scipy.sparse.matrix_csr, etc. + Matrix in the type of the underlying libraries. + """ + return _data.extract(self._data, format, copy) + + def contract(self, inplace: bool = False) -> Qobj: + """ + Contract subspaces of the tensor structure which are 1D. Not defined + on superoperators. If all dimensions are scalar, a Qobj of dimension + [[1], [1]] is returned, i.e. _multiple_ scalar dimensions are + contracted, but one is left. + + Parameters + ---------- + inplace: bool, optional + If ``True``, modify the dimensions in place. If ``False``, return + a copied object. + + Returns + ------- + out: :class:`.Qobj` + Quantum object with dimensions contracted. Will be ``self`` if + ``inplace`` is ``True``. + """ + if self.isket: + sub = [x for x in self.dims[0] if x > 1] or [1] + dims = [sub, [1]*len(sub)] + elif self.isbra: + sub = [x for x in self.dims[1] if x > 1] or [1] + dims = [[1]*len(sub), sub] + elif self.isoper or self.isoperket or self.isoperbra: + if self.isoper: + oper_dims = self.dims + elif self.isoperket: + oper_dims = self.dims[0] + else: + oper_dims = self.dims[1] + if len(oper_dims[0]) != len(oper_dims[1]): + raise ValueError("cannot parse Qobj dimensions: " + + repr(self.dims)) + dims_ = [ + (x, y) for x, y in zip(oper_dims[0], oper_dims[1]) + if x > 1 or y > 1 + ] or [(1, 1)] + dims = [[x for x, _ in dims_], [y for _, y in dims_]] + if self.isoperket: + dims = [dims, [1]] + elif self.isoperbra: + dims = [[1], dims] + else: + raise TypeError("not defined for superoperators") + if inplace: + self.dims = dims + return self + return Qobj(self.data.copy(), dims=dims, copy=False) + + + def permute(self, order: list) -> Qobj: + """ + Permute the tensor structure of a quantum object. For example, + + ``qutip.tensor(x, y).permute([1, 0])`` + + will give the same result as + + ``qutip.tensor(y, x)`` + + and + + ``qutip.tensor(a, b, c).permute([1, 2, 0])`` + + will be the same as + + ``qutip.tensor(b, c, a)`` + + For regular objects (bras, kets and operators) we expect ``order`` to + be a flat list of integers, which specifies the new order of the tensor + product. + + For superoperators, we expect ``order`` to be something like + + ``[[0, 2], [1, 3]]`` + + which tells us to permute according to [0, 2, 1, 3], and then group + indices according to the length of each sublist. As another example, + permuting a superoperator with dimensions of + + ``[[[1, 2, 3], [1, 2, 3]], [[1, 2, 3], [1, 2, 3]]]`` + + by an ``order`` + + ``[[0, 3], [1, 4], [2, 5]]`` + + should give a new object with dimensions + + ``[[[1, 1], [2, 2], [3, 3]], [[1, 1], [2, 2], [3, 3]]]``. + + Parameters + ---------- + order : list + List of indices specifying the new tensor order. + + Returns + ------- + P : :class:`.Qobj` + Permuted quantum object. + """ + if self.type in ('bra', 'ket', 'oper'): + structure = self.dims[1] if self.isbra else self.dims[0] + new_structure = [structure[x] for x in order] + if self.isbra: + dims = [self.dims[0], new_structure] + elif self.isket: + dims = [new_structure, self.dims[1]] + else: + if self._dims[0] != self._dims[1]: + raise TypeError("undefined for non-square operators") + dims = [new_structure, new_structure] + data = _data.permute.dimensions(self.data, structure, order) + return Qobj(data, + dims=dims, + isherm=self._isherm, + isunitary=self._isunitary, + copy=False) + # If we've got here, we're some form of superoperator, so we work with + # the flattened structure. + flat_order = flatten(order) + flat_structure = flatten(self.dims[1] if self.isoperbra + else self.dims[0]) + new_structure = unflatten([flat_structure[x] for x in flat_order], + enumerate_flat(order)) + if self.isoperbra: + dims = [self.dims[0], new_structure] + elif self.isoperket: + dims = [new_structure, self.dims[1]] + else: + if self._dims[0] != self._dims[1]: + raise TypeError("undefined for non-square operators") + dims = [new_structure, new_structure] + data = _data.permute.dimensions(self.data, flat_structure, flat_order) + return Qobj(data, + dims=dims, + superrep=self.superrep, + copy=False) + + def tidyup(self, atol: float = None) -> Qobj: + """ + Removes small elements from the quantum object. + + Parameters + ---------- + atol : float + Absolute tolerance used by tidyup. Default is set + via qutip global settings parameters. + + Returns + ------- + oper : :class:`.Qobj` + Quantum object with small elements removed. + """ + atol = atol or settings.core['auto_tidyup_atol'] + self.data = _data.tidyup(self.data, atol) + return self + + def transform( + self, + inpt: list[Qobj] | ArrayLike, + inverse: bool = False + ) -> Qobj: + """Basis transform defined by input array. + + Input array can be a ``matrix`` defining the transformation, + or a ``list`` of kets that defines the new basis. + + Parameters + ---------- + inpt : array_like + A ``matrix`` or ``list`` of kets defining the transformation. + inverse : bool + Whether to return inverse transformation. + + Returns + ------- + oper : :class:`.Qobj` + Operator in new basis. + + Notes + ----- + This function is still in development. + """ + # TODO: What about superoper + # TODO: split into each cases. + if isinstance(inpt, list) or (isinstance(inpt, np.ndarray) and + inpt.ndim == 1): + if len(inpt) != max(self.shape): + raise TypeError( + 'Invalid size of ket list for basis transformation') + base = np.hstack([psi.full() for psi in inpt]) + S = _data.adjoint(_data.create(base)) + elif isinstance(inpt, Qobj) and inpt.isoper: + S = inpt.data + elif isinstance(inpt, np.ndarray): + S = _data.create(inpt).conj() + else: + raise TypeError('Invalid operand for basis transformation') + + # transform data + if inverse: + if self.isket: + data = _data.matmul(S.adjoint(), self.data) + elif self.isbra: + data = _data.matmul(self.data, S) + else: + data = _data.matmul(_data.matmul(S.adjoint(), self.data), S) + else: + if self.isket: + data = _data.matmul(S, self.data) + elif self.isbra: + data = _data.matmul(self.data, S.adjoint()) + else: + data = _data.matmul(_data.matmul(S, self.data), S.adjoint()) + return Qobj(data, + dims=self.dims, + isherm=self._isherm, + superrep=self.superrep, + copy=False) + + # TODO: split + def overlap(self, other: Qobj) -> complex: + """ + Overlap between two state vectors or two operators. + + Gives the overlap (inner product) between the current bra or ket Qobj + and and another bra or ket Qobj. It gives the Hilbert-Schmidt overlap + when one of the Qobj is an operator/density matrix. + + Parameters + ---------- + other : :class:`.Qobj` + Quantum object for a state vector of type 'ket', 'bra' or density + matrix. + + Returns + ------- + overlap : complex + Complex valued overlap. + + Raises + ------ + TypeError + Can only calculate overlap between a bra, ket and density matrix + quantum objects. + """ + if not isinstance(other, Qobj): + raise TypeError("".join([ + "cannot calculate overlap with non-quantum object ", + repr(other), + ])) + if ( + self.type not in ('ket', 'bra', 'oper') + or other.type not in ('ket', 'bra', 'oper') + ): + msg = "only bras, kets and density matrices have defined overlaps" + raise TypeError(msg) + left, right = self._data.adjoint(), other.data + if self.isoper or other.isoper: + if not self.isoper: + left = _data.project(left) + if not other.isoper: + right = _data.project(right) + return _data.trace(_data.matmul(left, right)) + if other.isbra: + right = right.adjoint() + out = _data.inner(left, right, self.isket) + if self.isket and other.isbra: + # In this particular case, we've basically doing + # conj(other.overlap(self)) + # so we take care to conjugate the output. + out = np.conj(out) + return out + + @property + def ishp(self) -> bool: + return False + + @property + def iscp(self) -> bool: + return False + + @property + def istp(self) -> bool: + return False + + @property + def iscptp(self) -> bool: + return False + + @property + def isherm(self) -> bool: + return False + + @property + def _isherm(self) -> bool: + # Weak version of `isherm`, does not compute if missing + return self._flags.get("isherm", None) + + @property + def isunitary(self) -> bool: + return False + + @property + def _isunitary(self) -> bool: + # Weak version of `isunitary`, does not compute if missing + return self._flags.get("isunitary", None) + + @property + def shape(self) -> tuple[int, int]: + """Return the shape of the Qobj data.""" + return self._data.shape + + @property + def isoper(self) -> bool: + """Indicates if the Qobj represents an operator.""" + return False + + @property + def isbra(self) -> bool: + """Indicates if the Qobj represents a bra state.""" + return False + + @property + def isket(self) -> bool: + """Indicates if the Qobj represents a ket state.""" + return False + + @property + def issuper(self) -> bool: + """Indicates if the Qobj represents a superoperator.""" + return False + + @property + def isoperket(self) -> bool: + """Indicates if the Qobj represents a operator-ket state.""" + return False + + @property + def isoperbra(self) -> bool: + """Indicates if the Qobj represents a operator-bra state.""" + return False + + @property + def superrep(self) -> str: + return self._dims.superrep + + @superrep.setter + def superrep(self, super_rep: str): + self._dims = self._dims.replace_superrep(super_rep) + + def ptrace(self, sel: int | list[int], dtype: LayerType = None) -> Qobj: + """ + Take the partial trace of the quantum object leaving the selected + subspaces. In other words, trace out all subspaces which are _not_ + passed. + + This is typically a function which acts on operators; bras and kets + will be promoted to density matrices before the operation takes place + since the partial trace is inherently undefined on pure states. + + For operators which are currently being represented as states in the + superoperator formalism (i.e. the object has type `operator-ket` or + `operator-bra`), the partial trace is applied as if the operator were + in the conventional form. This means that for any operator `x`, + ``operator_to_vector(x).ptrace(0) == operator_to_vector(x.ptrace(0))`` + and similar for `operator-bra`. + + The story is different for full superoperators. In the formalism that + QuTiP uses, if an operator has dimensions (`dims`) of + `[[2, 3], [2, 3]]` then it can be represented as a state on a Hilbert + space of dimensions `[2, 3, 2, 3]`, and a superoperator would be an + operator which acts on this joint space. This function performs the + partial trace on superoperators by letting the selected components + refer to elements of the _joint_ _space_, and then returns a regular + operator (of type `oper`). + + Parameters + ---------- + sel : int or iterable of int + An ``int`` or ``list`` of components to keep after partial trace. + The selected subspaces will _not_ be reordered, no matter order + they are supplied to `ptrace`. + + dtype : type, str + The matrix format of output. + + Returns + ------- + oper : :class:`.Qobj` + Quantum object representing partial trace with selected components + remaining. + """ + return qutip.ptrace(self, sel, dtype) + + def purity(self) -> complex: + """Calculate purity of a quantum object. + + Returns + ------- + state_purity : float + Returns the purity of a quantum object. + For a pure state, the purity is 1. + For a mixed state of dimension `d`, 1/d<=purity<1. + + """ + if self.type in ("super", "operator-ket", "operator-bra"): + raise TypeError('purity is only defined for states.') + if self.isket or self.isbra: + return _data.norm.l2(self._data)**2 + return _data.trace(_data.matmul(self._data, self._data)).real diff --git a/qutip/core/qobj/_operator.py b/qutip/core/qobj/_operator.py new file mode 100644 index 0000000000..25d5a67060 --- /dev/null +++ b/qutip/core/qobj/_operator.py @@ -0,0 +1,711 @@ +from ._base import Qobj, _QobjBuilder +from ...settings import settings +import numpy as np +from qutip.typing import LayerType +import qutip +from typing import Any, Literal +import numbers +from .. import data as _data +from ..dimensions import enumerate_flat, collapse_dims_super +import warnings +from typing import Any, Literal + + +__all__ = [] + +class Operator(Qobj): + def __init__(self, data, dims, **flags): + super().__init__(data, dims, **flags) + if self._dims.type not in ["oper"]: + raise ValueError( + f"Expected operator dimensions, but got {self._dims.type}" + ) + + @property + def isoper(self) -> bool: + return True + + @property + def isherm(self) -> bool: + if self._flags.get("isherm", None) is None: + self._flags["isherm"] = _data.isherm(self._data) + return self._flags["isherm"] + + @isherm.setter + def isherm(self, isherm: bool): + self._flags["isherm"] = isherm + + @property + def isunitary(self) -> bool: + if self._flags.get("isunitary", None) is None: + if not self.isoper or self._data.shape[0] != self._data.shape[1]: + self._flags["isunitary"] = False + else: + cmp = _data.matmul(self._data, self._data.adjoint()) + iden = _data.identity_like(cmp) + self._flags["isunitary"] = _data.iszero( + _data.sub(cmp, iden), tol=settings.core['atol'] + ) + return self._flags["isunitary"] + + @isunitary.setter + def isunitary(self, isunitary: bool): + self._flags["isunitary"] = isunitary + + @property + def ishp(self) -> bool: + return True + + @property + def iscp(self) -> bool: + return True + + @property + def istp(self) -> bool: + return self.isunitary + + @property + def iscptp(self) -> bool: + return self.isunitary + + def matrix_element(self, bra: Qobj, ket: Qobj) -> Qobj: + """Calculates a matrix element. + + Gives the matrix element for the quantum object sandwiched between a + `bra` and `ket` vector. + + Parameters + ---------- + bra : :class:`.Qobj` + Quantum object of type 'bra' or 'ket' + + ket : :class:`.Qobj` + Quantum object of type 'ket'. + + Returns + ------- + elem : complex + Complex valued matrix element. + + Notes + ----- + It is slightly more computationally efficient to use a ket + vector for the 'bra' input. + + """ + if bra.type not in ('bra', 'ket') or ket.type not in ('bra', 'ket'): + msg = "Can only calculate matrix elements between a bra and a ket." + raise TypeError(msg) + left, op, right = bra.data, self.data, ket.data + if ket.isbra: + right = right.adjoint() + return _data.inner_op(left, op, right, bra.isket) + + def dnorm(self, B: Qobj = None) -> float: + """Calculates the diamond norm, or the diamond distance to another + operator. + + Parameters + ---------- + B : :class:`.Qobj` or None + If B is not None, the diamond distance d(A, B) = dnorm(A - B) + between this operator and B is returned instead of the diamond norm. + + Returns + ------- + d : float + Either the diamond norm of this operator, or the diamond distance + from this operator to B. + + """ + return qutip.dnorm(self, B) + + def __call__(self, other: Qobj) -> Qobj: + """ + Acts this Qobj on another Qobj either by left-multiplication, + or by vectorization and devectorization, as + appropriate. + """ + if not isinstance(other, Qobj): + raise TypeError("Only defined for quantum objects.") + if other.type not in ["ket"]: + raise TypeError(self.type + " cannot act on " + other.type) + return self.__matmul__(other) + + def __pow__(self, n: int, m=None) -> Qobj: + # calculates powers of Qobj + if ( + m is not None + or not isinstance(n, numbers.Integral) + or n < 0 + ): + return NotImplemented + return Qobj(_data.pow(self._data, n), + dims=self._dims, + isherm=self._isherm, + isunitary=self._isunitary, + copy=False) + + def tr(self) -> complex: + """Trace of a quantum object. + + Returns + ------- + trace : float + Returns the trace of the quantum object. + + """ + out = _data.trace(self._data) + # This ensures that trace can return something that is not a number such + # as a `tensorflow.Tensor` in qutip-tensorflow. + if settings.core["auto_real_casting"] and self.isherm: + out = out.real + return out + + def diag(self) -> np.ndarray: + """Diagonal elements of quantum object. + + Returns + ------- + diags : array + Returns array of ``real`` values if operators is Hermitian, + otherwise ``complex`` values are returned. + """ + # TODO: add a `diagonal` method to the data layer? + out = _data.to(_data.CSR, self.data).as_scipy().diagonal() + if settings.core["auto_real_casting"] and self.isherm: + out = np.real(out) + return out + + def expm(self, dtype: LayerType = None) -> Qobj: + """Matrix exponential of quantum operator. + + Input operator must be square. + + Parameters + ---------- + dtype : type + The data-layer type that should be output. + + Returns + ------- + oper : :class:`.Qobj` + Exponentiated quantum operator. + + Raises + ------ + TypeError + Quantum operator is not square. + """ + if self.shape[0] != self.shape[1]: + raise TypeError("expm is only valid for square operators") + if dtype is None and isinstance(self.data, (_data.CSR, _data.Dia)): + dtype = _data.Dense + return Qobj(_data.expm(self._data, dtype=dtype), + dims=self._dims, + isherm=self._isherm, + copy=False) + + def logm(self) -> Qobj: + """Matrix logarithm of quantum operator. + + Input operator must be square. + + Returns + ------- + oper : :class:`.Qobj` + Logarithm of the quantum operator. + + Raises + ------ + TypeError + Quantum operator is not square. + """ + if self.shape[0] != self.shape[1]: + raise TypeError("logm is only valid for square operators") + return Qobj(_data.logm(self._data), + dims=self._dims, + isherm=self._isherm, + copy=False) + + def sqrtm( + self, + sparse: bool = False, + tol: float = 0, + maxiter: int = 100000 + ) -> Qobj: + """ + Sqrt of a quantum operator. Operator must be square. + + Parameters + ---------- + sparse : bool + Use sparse eigenvalue/vector solver. + tol : float + Tolerance used by sparse solver (0 = machine precision). + maxiter : int + Maximum number of iterations used by sparse solver. + + Returns + ------- + oper : :class:`.Qobj` + Matrix square root of operator. + + Raises + ------ + TypeError + Quantum object is not square. + + Notes + ----- + The sparse eigensolver is much slower than the dense version. + Use sparse only if memory requirements demand it. + """ + if self.shape[0] != self.shape[1]: + raise TypeError('sqrt only valid on square matrices') + return Qobj(_data.sqrtm(self._data), + dims=self._dims, + copy=False) + + def cosm(self) -> Qobj: + """Cosine of a quantum operator. + + Operator must be square. + + Returns + ------- + oper : :class:`.Qobj` + Matrix cosine of operator. + + Raises + ------ + TypeError + Quantum object is not square. + + Notes + ----- + Uses the Q.expm() method. + + """ + if self.shape[0] != self.shape[1]: + raise TypeError('invalid operand for matrix cosine') + return 0.5 * ((1j * self).expm() + (-1j * self).expm()) + + def sinm(self) -> Qobj: + """Sine of a quantum operator. + + Operator must be square. + + Returns + ------- + oper : :class:`.Qobj` + Matrix sine of operator. + + Raises + ------ + TypeError + Quantum object is not square. + + Notes + ----- + Uses the Q.expm() method. + """ + if self.shape[0] != self.shape[1]: + raise TypeError('invalid operand for matrix sine') + return -0.5j * ((1j * self).expm() - (-1j * self).expm()) + + def inv(self, sparse: bool = False) -> Qobj: + """Matrix inverse of a quantum operator + + Operator must be square. + + Returns + ------- + oper : :class:`.Qobj` + Matrix inverse of operator. + + Raises + ------ + TypeError + Quantum object is not square. + """ + if self.shape[0] != self.shape[1]: + raise TypeError('Invalid operand for matrix inverse') + if isinstance(self.data, _data.CSR) and not sparse: + data = _data.to(_data.Dense, self.data) + else: + data = self.data + + return Qobj(_data.inv(data), + dims=[self._dims[1], self._dims[0]], + copy=False) + + def check_herm(self) -> bool: + """Check if the quantum object is hermitian. + + Returns + ------- + isherm : bool + Returns the new value of isherm property. + """ + self.isherm = None + return self.isherm + + def eigenstates( + self, + sparse: bool = False, + sort: Literal["low", "high"] = 'low', + eigvals: int = 0, + tol: float = 0, + maxiter: int = 100000, + phase_fix: int = None + ) -> tuple[np.ndarray, list[Qobj]]: + """Eigenstates and eigenenergies. + + Eigenstates and eigenenergies are defined for operators and + superoperators only. + + Parameters + ---------- + sparse : bool + Use sparse Eigensolver + + sort : str + Sort eigenvalues (and vectors) 'low' to high, or 'high' to low. + + eigvals : int + Number of requested eigenvalues. Default is all eigenvalues. + + tol : float + Tolerance used by sparse Eigensolver (0 = machine precision). + The sparse solver may not converge if the tolerance is set too low. + + maxiter : int + Maximum number of iterations performed by sparse solver (if used). + + phase_fix : int, None + If not None, set the phase of each kets so that ket[phase_fix,0] + is real positive. + + Returns + ------- + eigvals : array + Array of eigenvalues for operator. + + eigvecs : array + Array of quantum operators representing the oprator eigenkets. + Order of eigenkets is determined by order of eigenvalues. + + Notes + ----- + The sparse eigensolver is much slower than the dense version. + Use sparse only if memory requirements demand it. + + """ + if self.shape[0] != self.shape[1]: + raise TypeError('eigenstates only valid on square matrices') + if isinstance(self.data, _data.CSR) and sparse: + evals, evecs = _data.eigs_csr(self.data, + isherm=self._isherm, + sort=sort, eigvals=eigvals, tol=tol, + maxiter=maxiter) + elif isinstance(self.data, (_data.CSR, _data.Dia)): + evals, evecs = _data.eigs(_data.to(_data.Dense, self.data), + isherm=self._isherm, + sort=sort, eigvals=eigvals) + else: + evals, evecs = _data.eigs(self.data, isherm=self._isherm, + sort=sort, eigvals=eigvals) + + if self.type == 'super': + new_dims = [self.dims[0], [1]] + else: + new_dims = [self.dims[0], [1]*len(self.dims[0])] + ekets = np.empty((evecs.shape[1],), dtype=object) + ekets[:] = [Qobj(vec, dims=new_dims, copy=False) + for vec in _data.split_columns(evecs, False)] + norms = np.array([ket.norm() for ket in ekets]) + if phase_fix is None: + phase = np.array([1] * len(ekets)) + else: + phase = np.array([np.abs(ket[phase_fix, 0]) / ket[phase_fix, 0] + if ket[phase_fix, 0] else 1 + for ket in ekets]) + return evals, ekets / norms * phase + + def eigenenergies( + self, + sparse: bool = False, + sort: Literal["low", "high"] = 'low', + eigvals: int = 0, + tol: float = 0, + maxiter: int = 100000, + ) -> np.ndarray: + """Eigenenergies of a quantum object. + + Eigenenergies (eigenvalues) are defined for operators or superoperators + only. + + Parameters + ---------- + sparse : bool + Use sparse Eigensolver + sort : str + Sort eigenvalues 'low' to high, or 'high' to low. + eigvals : int + Number of requested eigenvalues. Default is all eigenvalues. + tol : float + Tolerance used by sparse Eigensolver (0=machine precision). + The sparse solver may not converge if the tolerance is set too low. + maxiter : int + Maximum number of iterations performed by sparse solver (if used). + + Returns + ------- + eigvals : array + Array of eigenvalues for operator. + + Notes + ----- + The sparse eigensolver is much slower than the dense version. + Use sparse only if memory requirements demand it. + + """ + if self.shape[0] != self.shape[1]: + raise TypeError('eigenenergies only valid on square matrices') + # TODO: consider another way of handling the dispatch here. + if isinstance(self.data, _data.CSR) and sparse: + return _data.eigs_csr(self.data, + vecs=False, + isherm=self._isherm, + sort=sort, eigvals=eigvals, + tol=tol, maxiter=maxiter) + elif isinstance(self.data, (_data.CSR, _data.Dia)): + return _data.eigs(_data.to(_data.Dense, self.data), + vecs=False, isherm=self._isherm, + sort=sort, eigvals=eigvals) + + return _data.eigs(self.data, + vecs=False, + isherm=self._isherm, sort=sort, eigvals=eigvals) + + def groundstate( + self, + sparse: bool = False, + tol: float = 0, + maxiter: int = 100000, + safe: bool = True + ) -> tuple[float, Qobj]: + """Ground state Eigenvalue and Eigenvector. + + Defined for quantum operators only. + + Parameters + ---------- + sparse : bool + Use sparse Eigensolver + tol : float + Tolerance used by sparse Eigensolver (0 = machine precision). + The sparse solver may not converge if the tolerance is set too low. + maxiter : int + Maximum number of iterations performed by sparse solver (if used). + safe : bool (default=True) + Check for degenerate ground state + + Returns + ------- + eigval : float + Eigenvalue for the ground state of quantum operator. + eigvec : :class:`.Qobj` + Eigenket for the ground state of quantum operator. + + Notes + ----- + The sparse eigensolver is much slower than the dense version. + Use sparse only if memory requirements demand it. + """ + if self.shape[0] != self.shape[1]: + raise TypeError('groundstate only valid on square matrices') + eigvals = 2 if safe else 1 + evals, evecs = self.eigenstates(sparse=sparse, eigvals=eigvals, + tol=tol, maxiter=maxiter) + + if safe: + tol = tol or settings.core['atol'] + # This tol should be less strick than the tol for the eigensolver + # so it's numerical errors are not seens as degenerate states. + if (evals[1]-evals[0]) <= 10 * tol: + warnings.warn("Ground state may be degenerate.", UserWarning) + return evals[0], evecs[0] + + def trunc_neg(self, method: Literal["clip", "sgs"] = "clip") -> Qobj: + """Truncates negative eigenvalues and renormalizes. + + Returns a new Qobj by removing the negative eigenvalues + of this instance, then renormalizing to obtain a valid density + operator. + + Parameters + ---------- + method : str + Algorithm to use to remove negative eigenvalues. "clip" + simply discards negative eigenvalues, then renormalizes. + "sgs" uses the SGS algorithm (doi:10/bb76) to find the + positive operator that is nearest in the Shatten 2-norm. + + Returns + ------- + oper : :class:`.Qobj` + A valid density operator. + """ + if self.shape[0] != self.shape[1]: + raise TypeError('trunc_neg only valid on square matrices') + if not self.isherm: + raise ValueError("Must be a Hermitian operator to remove negative " + "eigenvalues.") + if method not in ('clip', 'sgs'): + raise ValueError("Method {} not recognized.".format(method)) + + eigvals, eigstates = self.eigenstates() + if all(eigval >= 0 for eigval in eigvals): + # All positive, so just renormalize. + return self.unit() + idx_nonzero = eigvals != 0 + eigvals = eigvals[idx_nonzero] + eigstates = eigstates[idx_nonzero] + + if method == 'clip': + eigvals[eigvals < 0] = 0 + elif method == 'sgs': + eigvals = eigvals[::-1] + eigstates = eigstates[::-1] + acc = 0.0 + n_eigs = len(eigvals) + for idx in reversed(range(n_eigs)): + if eigvals[idx] + acc / (idx + 1) >= 0: + break + acc += eigvals[idx] + eigvals[idx] = 0.0 + eigvals[:idx+1] += acc / (idx + 1) + out_data = _data.zeros(*self.shape) + for value, state in zip(eigvals, eigstates): + if value: + # add in 3-argument form is fused-add-multiply + out_data = _data.add(out_data, + _data.project(state.data), + value) + out_data = _data.mul(out_data, 1/_data.norm.trace(out_data)) + return Qobj(out_data, dims=self._dims, isherm=True, copy=False) + + def norm( + self, + norm: Literal["max", "fro", "tr", "one"] = "tr", + kwargs: dict[str, Any] = None + ) -> float: + """ + Norm of a quantum object. + + Default norm is the trace-norm. Other operator norms may be + specified using the `norm` parameter. + + Parameters + ---------- + norm : str, default: "tr" + Which type of norm to use. Allowed values are 'tr' for the trace + norm, 'fro' for the Frobenius norm, 'one' and 'max'. + + kwargs : dict, optional + Additional keyword arguments to pass on to the relevant norm + solver. See details for each norm function in :mod:`.data.norm`. + + Returns + ------- + norm : float + The requested norm of the operator or state quantum object. + """ + norm = norm or "tr" + if norm not in {'tr', 'fro', 'one', 'max'}: + raise ValueError( + "matrix norm must be in {'tr', 'fro', 'one', 'max'}" + ) + + kwargs = kwargs or {} + return { + 'tr': _data.norm.trace, + 'one': _data.norm.one, + 'max': _data.norm.max, + 'fro': _data.norm.frobenius, + }[norm](self._data, **kwargs) + + def unit( + self, + inplace: bool = False, + norm: Literal["l2", "max", "fro", "tr", "one"] = "tr", + kwargs: dict[str, Any] = None + ) -> Qobj: + """ + Operator or state normalized to unity. Uses norm from Qobj.norm(). + + Parameters + ---------- + inplace : bool + Do an in-place normalization + norm : str + Requested norm for states / operators. + kwargs : dict + Additional key-word arguments to be passed on to the relevant norm + function (see :meth:`.norm` for more details). + + Returns + ------- + obj : :class:`.Qobj` + Normalized quantum object. Will be the `self` object if in place. + """ + norm_ = self.norm(norm=norm, kwargs=kwargs) + if inplace: + self.data = _data.mul(self.data, 1 / norm_) + self.isherm = self._isherm if norm_.imag == 0 else None + self.isunitary = ( + self._isunitary + if abs(norm_) - 1 < settings.core['atol'] + else None + ) + out = self + else: + out = self / norm_ + return out + + +class Scalar(Operator): + # Scalar can be anything + def __init__(self, data, dims, **flags): + Qobj.__init__(self, data, dims, **flags) + if self._dims.type not in ["scalar"]: + raise ValueError( + f"Expected scalar dimensions, but got {self._dims.type}" + ) + + @property + def issuper(self) -> bool: + """Indicates if the Qobj represents a superoperator.""" + return self._dims.issuper + + @property + def isket(self) -> bool: + return not self.issuper + + @property + def isbra(self) -> bool: + return not self.issuper + + @property + def isoperket(self) -> bool: + return self.issuper + + @property + def isoperbra(self) -> bool: + return self.issuper + + +_QobjBuilder.qobjtype_to_class["scalar"] = Scalar +_QobjBuilder.qobjtype_to_class["oper"] = Operator diff --git a/qutip/core/qobj/_state.py b/qutip/core/qobj/_state.py new file mode 100644 index 0000000000..b7527fcfba --- /dev/null +++ b/qutip/core/qobj/_state.py @@ -0,0 +1,153 @@ +from ._base import Qobj, _QobjBuilder +from .. import data as _data +from typing import Any, Literal + +__all__ = [] + +class _StateQobj(Qobj): + def proj(self) -> Qobj: + """Form the projector from a given ket or bra vector. + + Parameters + ---------- + Q : :class:`.Qobj` + Input bra or ket vector + + Returns + ------- + P : :class:`.Qobj` + Projection operator. + """ + dims = ([self._dims[0], self._dims[0]] if self.isket + else [self._dims[1], self._dims[1]]) + return Qobj(_data.project(self._data), + dims=dims, + isherm=True, + copy=False) + + def norm( + self, + norm: Literal["l2", "max"] = "l2", + kwargs: dict[str, Any] = None + ) -> float: + """ + Norm of a quantum object. + + Default norm is L2-norm. Other + ket and operator norms may be specified using the `norm` parameter. + + Parameters + ---------- + norm : str, default: "l2" + Which type of norm to use. Allowed values are 'l2' and 'max'. + + kwargs : dict, optional + Additional keyword arguments to pass on to the relevant norm + solver. See details for each norm function in :mod:`.data.norm`. + + Returns + ------- + norm : float + The requested norm of the operator or state quantum object. + """ + norm = norm or "l2" + if norm not in ["l2", "max"]: + raise ValueError( + "vector norm must be in {'k2', 'max'}" + ) + + kwargs = kwargs or {} + return { + 'l2': _data.norm.l2, + 'max': _data.norm.max + }[norm](self._data, **kwargs) + + def unit( + self, + inplace: bool = False, + norm: Literal["l2", "max"] = "l2", + kwargs: dict[str, Any] = None + ) -> Qobj: + """ + Operator or state normalized to unity. Uses norm from Qobj.norm(). + + Parameters + ---------- + inplace : bool, default: False + Do an in-place normalization + norm : str, default: "l2" + Requested norm for states / operators. + kwargs : dict, optional + Additional key-word arguments to be passed on to the relevant norm + function (see :meth:`.norm` for more details). + + Returns + ------- + obj : :class:`.Qobj` + Normalized quantum object. Will be the `self` object if in place. + """ + norm_ = self.norm(norm=norm, kwargs=kwargs) + if inplace: + self.data = _data.mul(self.data, 1 / norm_) + out = self + else: + out = self / norm_ + return out + + +class Bra(_StateQobj): + def __init__(self, data, dims, **flags): + super().__init__(data, dims, **flags) + if self._dims.type != "bra": + raise ValueError( + f"Expected bra dimensions, but got {self._dims.type}" + ) + + @property + def isbra(self) -> bool: + return True + + +class Ket(_StateQobj): + def __init__(self, data, dims, **flags): + super().__init__(data, dims, **flags) + if self._dims.type != "ket": + raise ValueError( + f"Expected ket dimensions, but got {self._dims.type}" + ) + + @property + def isket(self) -> bool: + return True + + +class OperKet(_StateQobj): + def __init__(self, data, dims, **flags): + super().__init__(data, dims, **flags) + if self._dims.type != "operator-ket": + raise ValueError( + f"Expected ket dimensions, but got {self._dims.type}" + ) + + @property + def isoperket(self) -> bool: + return True + + +class OperBra(_StateQobj): + def __init__(self, data, dims, **flags): + super().__init__(data, dims, **flags) + if self._dims.type != "operator-bra": + raise ValueError( + f"Expected ket dimensions, but got {self._dims.type}" + ) + + @property + def isoperbra(self) -> bool: + return True + + +_QobjBuilder.qobjtype_to_class["ket"] = Ket +_QobjBuilder.qobjtype_to_class["bra"] = Bra +_QobjBuilder.qobjtype_to_class["operator-ket"] = OperKet +_QobjBuilder.qobjtype_to_class["operator-bra"] = OperBra diff --git a/qutip/core/qobj/_super.py b/qutip/core/qobj/_super.py new file mode 100644 index 0000000000..8312f12b51 --- /dev/null +++ b/qutip/core/qobj/_super.py @@ -0,0 +1,587 @@ +from ._base import Qobj, _QobjBuilder, _require_equal_type +from ._operator import Operator +import qutip +from ..dimensions import enumerate_flat, collapse_dims_super, Dimensions +import numpy as np +from ...settings import settings +from .. import data as _data +from typing import Sequence +from qutip.typing import LayerType, DimensionLike + + +__all__ = ["KrausMap"] + + +class SuperOperator(Qobj): + def __init__(self, data, dims, **flags): + super().__init__(data, dims, **flags) + if self._dims.type not in ["super"]: + raise ValueError( + f"Expected super operator dimensions, but got {self._dims.type}" + ) + + @property + def issuper(self) -> bool: + return True + + @property + def isherm(self) -> bool: + if self._flags.get("isherm", None) is None: + self._flags["isherm"] = _data.isherm(self._data) + return self._flags["isherm"] + + @isherm.setter + def isherm(self, isherm: bool): + self._flags["isherm"] = isherm + + @property + def isunitary(self) -> bool: + if self._flags.get("isunitary", None) is None: + if not self.isoper or self._data.shape[0] != self._data.shape[1]: + self._flags["isunitary"] = False + else: + cmp = _data.matmul(self._data, self._data.adjoint()) + iden = _data.identity_like(cmp) + self._flags["isunitary"] = _data.iszero( + _data.sub(cmp, iden), tol=settings.core['atol'] + ) + return self._flags["isunitary"] + + @isunitary.setter + def isunitary(self, isunitary: bool): + self._flags["isunitary"] = isunitary + + @property + def ishp(self) -> bool: + if self._flags.get("ishp", None) is None: + try: + J = qutip.to_choi(self) + self._flags["ishp"] = J.isherm + except: + self._flags["ishp"] = False + + return self._flags["ishp"] + + @property + def iscp(self) -> bool: + if self._flags.get("iscp", None) is None: + # We can test with either Choi or chi, since the basis + # transformation between them is unitary and hence preserves + # the CP and TP conditions. + J = self if self.superrep in ('choi', 'chi') else qutip.to_choi(self) + # If J isn't hermitian, then that could indicate either that J is not + # normal, or is normal, but has complex eigenvalues. In either case, + # it makes no sense to then demand that the eigenvalues be + # non-negative. + self._flags["iscp"] = ( + J.isherm + and np.all(J.eigenenergies() >= -settings.core['atol']) + ) + return self._flags["iscp"] + + @property + def istp(self) -> bool: + if self._flags.get("istp", None) is None: + # Normalize to a super of type choi or chi. + # We can test with either Choi or chi, since the basis + # transformation between them is unitary and hence + # preserves the CP and TP conditions. + if self.issuper and self.superrep in ('choi', 'chi'): + qobj = self + else: + qobj = qutip.to_choi(self) + # Possibly collapse dims. + if any([len(index) > 1 + for super_index in qobj.dims + for index in super_index]): + qobj = Qobj(qobj.data, + dims=collapse_dims_super(qobj.dims), + superrep=qobj.superrep, + copy=False) + # We use the condition from John Watrous' lecture notes, + # Tr_1(J(Phi)) = identity_2. + # See: https://cs.uwaterloo.ca/~watrous/LectureNotes.html, + # Theory of Quantum Information (Fall 2011), theorem 5.4. + tr_oper = qobj.ptrace([0]) + self._flags["istp"] = np.allclose( + tr_oper.full(), + np.eye(tr_oper.shape[0]), + atol=settings.core['atol'] + ) + return self._flags["istp"] + + @property + def iscptp(self) -> bool: + if ( + self._flags.get("iscp", None) is None + and self._flags.get("istp", None) is None + ): + reps = ('choi', 'chi') + q_oper = qutip.to_choi(self) if self.superrep not in reps else self + self._flags["iscp"] = q_oper.iscp + self._flags["istp"] = q_oper.istp + return self.iscp and self.istp + + def dnorm(self, B: Qobj = None) -> float: + """Calculates the diamond norm, or the diamond distance to another + operator. + + Parameters + ---------- + B : :class:`.Qobj` or None + If B is not None, the diamond distance d(A, B) = dnorm(A - B) + between this operator and B is returned instead of the diamond norm. + + Returns + ------- + d : float + Either the diamond norm of this operator, or the diamond distance + from this operator to B. + + """ + return qutip.dnorm(self, B) + + def __call__(self, other: Qobj) -> Qobj: + """ + Acts this Qobj on another Qobj either by left-multiplication, + or by vectorization and devectorization, as + appropriate. + """ + if not isinstance(other, Qobj): + raise TypeError("Only defined for quantum objects.") + if other.type not in ["ket", "oper"]: + raise TypeError(self.type + " cannot act on " + other.type) + if other.isket: + other = other.proj() + return qutip.vector_to_operator(self @ qutip.operator_to_vector(other)) + + # Method from operator are merged on case per case basis + + def _warn(f, name=None): + import warnings + def func(*a, **kw): + warnings.warn(f"SuperOperator.{name} used") + return f(*a, **kw) + return func + + # Matrix operations, should we support them? + # Can't be easily applied on kraus map + # __pow__ = Operator.__pow__ + # expm = Operator.expm + # logm = Operator.logm + # cosm = Operator.cosm + # cosm = Operator.cosm + sqrtm = _warn(Operator.sqrtm, "sqrtm") # Fidelity used for choi + # inv = Operator.inv + + # These depend on the matrix representation. + tr = Operator.tr + # diag = Operator.diag + + eigenstates = Operator.eigenstates # Could be modified to return dm + eigenenergies = Operator.eigenenergies + # groundstate = Operator.groundstate # Useful for super operator? + + def dual_chan(self) -> Qobj: + """Dual channel of quantum object representing a completely positive + map. + """ + # Uses the technique of Johnston and Kribs (arXiv:1102.0948), which + # is only valid for completely positive maps. + if not self.iscp: + raise ValueError("Dual channels are only implemented for CP maps.") + J = qutip.to_choi(self) + tensor_idxs = enumerate_flat(J.dims) + J_dual = qutip.tensor_swap( + J, + *( + list(zip(tensor_idxs[0][1], tensor_idxs[0][0])) + + list(zip(tensor_idxs[1][1], tensor_idxs[1][0])) + ) + ).trans() + J_dual.superrep = 'choi' + return J_dual + + def to_choi(self): + return qutip.to_choi(self) + + def to_super(self): + return qutip.to_super(self) + + def to_kraus(self): + if self.ishp: + return qutip.to_kraus(self) + else: + from qutip.core.superop_reps import _generalized_kraus + pp = [(u, v.dag()) for u, v in _generalized_kraus(self.to_choi())] + return KrausMap.generalizedKraus(general_terms=pp) + + +class KrausMap(SuperOperator): + """ + Parameter + --------- + arg : list of Qobj + + """ + @staticmethod + def _check_oper(oper, dims): + if not isinstance(oper, Qobj) and oper.isoper: + raise TypeError("KrausMap componants must be operators") + if oper._dims != dims: + raise TypeError( + "KrausMap operators must all have the same dimensions" + ) + + def __init__(self, arg, dims=None, copy=False, **_): + self._data = None + self._kmap = [] + self._hpmap = [] + self._genmap = [] + dims_ = None + for oper in arg: + if dims_ is None: + dims_ = oper._dims + KrausMap._check_oper(oper, dims_) + self._kmap.append(oper.copy() if copy else oper) + + if dims_ is not None: + self._dims = Dimensions.from_prepost(dims_, dims_.transpose()) + if dims is not None and Dimensions(dims) != self._dims: + raise ValueError( + "Provided dimensions do not match operators" + ) + elif dims is not None: + # Dimensions but no operators, this result in a map that always + # a zeros state. + self._dims = Dimensions(dims) + else: + raise ValueError("Missing information to initialise KrausMap") + + self._flags = { + "ishp": True, + "isherm": False, + "isunitary": False, + } + self.superrep = "kraus" + + @classmethod + def generalizedKraus( + cls, + kraus_terms: Sequence[Operator] = (), + hp_terms: Sequence[Operator] = (), + general_terms: Sequence[tuple[Operator, Operator]] = (), + copy: bool = False, + ) -> "KrausMap": + """ + Create a generalized Kraus Map. + + G(rho) -> sum_i A[i] @ rho @ B[i] + + Internaly hermiticity preserving terms are kept appart: + + G(rho) -> + sum_i K[i] @ rho @ K[i].dag() + + sum_i (H[i] @ rho + rho @ H[i].dag()) + + sum_i A[i] @ rho @ B[i] + + Parameters + ---------- + kraus_terms: list of Qobj + Terms in the form: K[i] @ rho @ K[i].dag() + + hp_terms: list of Qobj + Terms in the form: H[i] @ rho + rho @ H[i].dag() + + general_terms: list of tuple of Qobj, Qobj + Terms in the form: A[i] @ rho @ B[i] + """ + out = KrausMap.__new__(KrausMap) + dims_pre = None + dims_post = None + self._data = None + self._kmap = [] + self._hpmap = [] + self._genmap = [] + + for oper in kraus_terms: + if dims_pre is None: + dims_pre = oper._dims + dims_post = oper._dims.tr() + KrausMap._check_oper(oper, dims_pre) + self._kmap.append(oper.copy() if copy else oper) + + for oper in hp_terms: + if dims_pre is None: + dims_pre = oper._dims + dims_post = oper._dims.tr() + KrausMap._check_oper(oper, dims_pre) + out._hpmap.append(oper.copy() if copy else oper) + + for pre, post in general_terms: + if dims_pre is None: + dims_pre = pre._dims + dims_post = post._dims + KrausMap._check_oper(pre, dims_pre) + KrausMap._check_oper(post, dims_post) + out._genmap.append( + (pre.copy(), post.copy()) + if copy else (pre, post) + ) + + if dims_pre is None: + raise ValueError( + "Can't initialise general kraus map from no operators" + ) + + self._dims = Dimensions.from_prepost(dims_pre, dims_post) + self._flags = { + "isherm": False, + "isunitary": False, + } + self.superrep = "kraus" + + def _pre_dims(self): + return Dimensions(self._dims[1].oper[0], self._dims[0].oper[0]) + + def _post_dims(self): + return Dimensions(self._dims[0].oper[1], self._dims[1].oper[1]) + + def __len__(self): + return len(self._kmap) + len(self._hpmap) + len(self._genmap) + + def __call__(self, other: Qobj) -> Qobj: + if not isinstance(other, Qobj): + raise TypeError("Only defined for quantum objects.") + if (self.type, other.type) not in _CALL_ALLOWED: + raise TypeError(self.type + " cannot act on " + other.type) + if other.isket: + other = other.proj() + out = qzeros(self._dims[1]) + for oper in self._kmap: + out += oper @ other @ oper.dag() + for oper in self._hpmap: + out += oper @ other + other @ oper.dag() + for pre, post in self._genmap: + out += pre @ other @ post + return out + + @property + def data(self): + raise ValueError("data not defined") + + @property + def shape(self): + return self._pre_dims().shape + + def dtype(self): + raise NotImplementedError + + def to(self, data_type: LayerType, copy: bool=False) -> Qobj: + raise NotImplementedError + + @_require_equal_type + def __add__(self, other): + if other == 0: + return self.copy() + return KrausMap.generalizedKraus( + kraus_terms = self._kmap + other._kmap, + hp_terms = self._hpmap + other._hpmap, + general_terms = self._genmap + other._genmap, + ) + + @_require_equal_type + def __sub__(self, other): + if other == 0: + return self.copy() + return self + other * -1 + + def __neg__(self): + return other * -1 + + def __mul__(self, other): + if isinstance(other, Qobj): + return self.__matmul__(other) + + if other == 0: + return KrausMap([], dims=self._dims) + + try: + if other.real == other and other >= 0: + sq = other.real ** 0.5 + kmap = [u * sq for u in self._kmap] + genmap = [] + else: + kmap = [] + genmap = [(u * other, u.dag()) for u in self._kmap] + if other.real == other: + hpmap = [u * other for u in self._hpmap] + else: + hpmap = [] + genmap += [ + (u * other, qutip.qeye_like(u)) + for u in self._hpmap + ] + genmap += [ + (qutip.qeye_like(u), u.dag() * other) + for u in self._hpmap + ] + genmap += [ + (u, v * other) + for u, v in self._genmap + ] + except TypeError: + return NotImplemented + + return KrausMap.generalizedKraus( + kraus_terms = kmap, + hp_terms = hpmap, + general_terms = genmap, + ) + + def __matmul__(self, other): + return NotImplemented + + def __eq__(self, other): + if self is other: + return True + if not isinstance(other, Qobj) or self._dims != other._dims: + return False + return NotImplementedError + + def __and__(self, other): + raise NotImplementedError + + def dag(self): + raise NotImplementedError + + def conj(self): + raise NotImplementedError + + def trans(self): + raise NotImplementedError + + def copy(self): + return KrausMap.generalizedKraus( + kraus_terms=self._kmap, + hp_terms=self._hpmap, + general_terms=self._genmap, + copy=True + ) + + def __str__(self): + out = self._str_header() + if self._genmap or self._hpmap: + type = "Generalized " + out + return out + + def _repr_latex_(self): + return self.__str__() + + def __getitem__(self, ind): + raise NotImplementedError + + def full(self): + shape = self.shape + if not (self._genmap or self._hpmap): + out = np.zeros((len(self._kmap), *shape), dtype=complex) + for i, oper in enumerate(self._kmap): + out[i] = oper.full() + else: + nk = len(self._kmap) + nhp = len(self._hpmap) + ng = len(self._genmap) + out = np.zeros((nk + 2 * nhp + ng, 2, *shape), dtype=complex) + for i, oper in enumerate(self._kmap): + out[i, 0] = oper.full() + out[i, 1] = oper.dag().full() + for i, oper in enumerate(self._hpmap): + out[2 * i + nk, 0] = oper.full() + out[2 * i + nk, 1] = qutip.qeye(self._post_dims()) + out[2 * i + nk, 0] = qutip.qeye(self._pre_dims()) + out[2 * i + nk, 1] = oper.dag().full() + for i, (u, v) in enumerate(self._genmap): + out[i + nk + 2 * nhp, 0] = u.full() + out[i + nk + 2 * nhp, 1] = v.full() + return out + + def data_as(self): + raise NotImplementedError + + def contract(self): + raise NotImplementedError + + def permute(self): + raise NotImplementedError + + def tidyup(self): + raise NotImplementedError + + def norm(self): + raise NotImplementedError + + def unit(self): + raise NotImplementedError + + def transform(self): + raise NotImplementedError + + def overlap(self): + raise NotImplementedError + + def ptrace(self): + raise NotImplementedError + + def purity(self): + raise NotImplementedError + + def to_super(self): + out = 0 + for oper in self._kmap: + out += qutip.to_super(oper) + for oper in self._hpmap: + out += qutip.spre(oper) + out += qutip.spost(oper.dag()) + for pre, post in self._genmap: + out += qutip.sprepost(pre, post) + out.superrep = "super" + return out + + def to_choi(self): + out = 0 + print(len(self._kmap), len(self._hpmap), len(self._genmap)) + for oper in self._kmap: + vec = qutip.operator_to_vector(oper) + out += vec @ vec.dag() + for oper in self._hpmap: + ones = qeye_like(oper) + out += qutip.operator_to_vector(oper) @ qutip.operator_to_vector(ones).dag() + out += qutip.operator_to_vector(ones) @ qutip.operator_to_vector(oper).dag() + for pre, post in self._genmap: + out += qutip.operator_to_vector(pre) @ qutip.operator_to_vector(post).trans() + print(out) + out.superrep = "choi" + return out + + @property + def iscp(self) -> bool: + if self._flags.get("iscp", None) is None: + if not self._hpmap and not self._genmap: + self._flags["iscp"] = True + else: + choi = self.to_choi() + self._flags["iscp"] = choi.iscp + return self._flags["iscp"] + + @property + def ishp(self) -> bool: + if self._flags.get("ishp", None) is None: + if not self._genmap: + self._flags["ishp"] = True + else: + choi = self.to_choi() + self._flags["ishp"] = choi.ishp + return self._flags["ishp"] + + +_QobjBuilder.qobjtype_to_class["super"] = SuperOperator diff --git a/qutip/core/qobj.py b/qutip/core/qobj_backup.py similarity index 100% rename from qutip/core/qobj.py rename to qutip/core/qobj_backup.py diff --git a/qutip/core/states.py b/qutip/core/states.py index 9e78bce737..472c6d6f09 100644 --- a/qutip/core/states.py +++ b/qutip/core/states.py @@ -536,7 +536,7 @@ def thermal_dm( raise ValueError("The method option can only take " "values 'operator' or 'analytic'") out = qdiags(diags, 0, dims=[[N], [N]], shape=(N, N), dtype=dtype) - out._isherm = True + out.isherm = True return out diff --git a/qutip/core/superop_reps.py b/qutip/core/superop_reps.py index 029aeab2cc..04fc3dd514 100644 --- a/qutip/core/superop_reps.py +++ b/qutip/core/superop_reps.py @@ -22,7 +22,7 @@ from .superoperator import stack_columns, unstack_columns, sprepost from .tensor import tensor from .dimensions import flatten -from .qobj import Qobj +from .qobj import Qobj, KrausMap from .operators import identity, sigmax, sigmay, sigmaz from .states import basis @@ -130,15 +130,18 @@ def _choi_to_kraus(q_oper, tol=1e-9): vals, vecs = q_oper.eigenstates() dims = [q_oper.dims[0][1], q_oper.dims[0][0]] shape = (np.prod(q_oper.dims[0][1]), np.prod(q_oper.dims[0][0])) - return [Qobj(_data.mul(unstack_columns(vec.data, shape=shape), np.sqrt(val)), - dims=dims, copy='False') - for val, vec in zip(vals, vecs) if abs(val) >= tol] + opers = [ + Qobj(_data.mul(unstack_columns(vec.data, shape=shape), np.sqrt(val)), + dims=dims, copy='False') + for val, vec in zip(vals, vecs) if abs(val) >= tol + ] + return KrausMap(opers) # Individual conversions from Kraus operators are public because the output # list of Kraus operators is not itself a quantum object. -def kraus_to_choi(kraus_ops: list[Qobj]) -> Qobj: +def kraus_to_choi(kraus_op: KrausMap) -> Qobj: r""" Convert a list of Kraus operators into Choi representation of the channel. @@ -149,8 +152,8 @@ def kraus_to_choi(kraus_ops: list[Qobj]) -> Qobj: Parameters ---------- - kraus_ops : list[Qobj] - The list of Kraus operators to be converted to Choi representation. + kraus_op : KrausMap + The kraus operators to be converted to Choi representation. Returns ------- @@ -158,23 +161,9 @@ def kraus_to_choi(kraus_ops: list[Qobj]) -> Qobj: A quantum object representing the same map as ``kraus_ops``, such that ``choi.superrep == "choi"``. """ - len_op = np.prod(kraus_ops[0].shape) - # If Kraus ops have dims [M, N] in qutip notation (act on [N, N] density - # matrix and produce [M, M] d.m.), Choi matrix Hilbert space will - # be [[M, N], [M, N]] because Choi Hilbert space - # is (output space) x (input space). - choi_dims = [kraus_ops[0].dims] * 2 - # transform a list of Qobj matrices list[sum_ij k_ij |i>> = sum_I k_I |I>> - kraus_vectors = np.asarray([ - np.reshape(kraus_op.full(), len_op, order="F") - for kraus_op in kraus_ops - ]) - # sum_{I} |k_I|^2 |I>>< Qobj: @@ -381,6 +370,8 @@ def to_choi(q_oper: Qobj) -> Qobj: return _super_tofrom_choi(q_oper) if q_oper.superrep == 'chi': return _chi_to_choi(q_oper) + elif q_oper.superrep == 'kraus': + return q_oper.to_choi() else: raise TypeError(q_oper.superrep) elif q_oper.type == 'oper': @@ -424,6 +415,8 @@ def to_chi(q_oper: Qobj) -> Qobj: return _choi_to_chi(q_oper) elif q_oper.superrep == 'super': return _choi_to_chi(to_choi(q_oper)) + elif q_oper.superrep == 'kraus': + return _choi_to_chi(q_oper.to_choi()) else: raise TypeError(q_oper.superrep) elif q_oper.type == 'oper': @@ -466,6 +459,8 @@ def to_super(q_oper: Qobj) -> Qobj: return _super_tofrom_choi(q_oper) elif q_oper.superrep == 'chi': return to_super(to_choi(q_oper)) + elif q_oper.superrep == 'kraus': + return q_oper.to_super() else: raise ValueError( "Unrecognized superrep '{}'.".format(q_oper.superrep)) @@ -506,11 +501,13 @@ def to_kraus(q_oper: Qobj, tol: float=1e-9) -> list[Qobj]: decomposed into Kraus operators. """ if q_oper.issuper: + if q_oper.superrep == "kraus": + return q_oper.copy() if q_oper.superrep != 'choi': q_oper = to_choi(q_oper) return _choi_to_kraus(q_oper, tol) elif q_oper.isoper: # Assume unitary - return [q_oper] + return KrausMap([q_oper]) raise TypeError( "Conversion of Qobj with type={0.type} " "and superrep={0.superrep} to Kraus decomposition not " diff --git a/qutip/core/superoperator.py b/qutip/core/superoperator.py index 706c2826b7..c072969c66 100644 --- a/qutip/core/superoperator.py +++ b/qutip/core/superoperator.py @@ -11,7 +11,7 @@ from .qobj import Qobj from .cy.qobjevo import QobjEvo from . import data as _data -from .dimensions import Compound, SuperSpace, Space +from .dimensions import Compound, SuperSpace, Space, Dimensions def _map_over_compound_operators(f): @@ -140,6 +140,18 @@ def liouvillian( copy=False) +def liouvillian_kraus(H: Qobj, c_ops: list[Qobj]): + """ + TODO + """ + out = KrausMap.generalizedKraus( + kraus_terms = c_ops, + hp_terms = [-1.0j * H] + [c.dag() @ c * -0.5 for c in c_ops], + ) + out._flags["istp"] = False + return out + + @overload def lindblad_dissipator( a: Qobj, @@ -222,6 +234,18 @@ def lindblad_dissipator( return D.data if data_only else D +def lindblad_dissipator_kraus(H: Qobj, c_ops: list[Qobj]): + """ + TODO + """ + out = KrausMap.generalizedKraus( + kraus_terms = c_ops, + hp_terms = [c.dag() @ c * -0.5 for c in c_ops], + ) + out._flags["istp"] = False + return out + + @_map_over_compound_operators def operator_to_vector(op: Qobj) -> Qobj: """ @@ -245,7 +269,7 @@ def operator_to_vector(op: Qobj) -> Qobj: raise TypeError("Cannot convert object already " "in super representation") return Qobj(stack_columns(op.data), - dims=[op.dims, [1]], + dims=[op._dims, [1]], superrep="super", copy=False) @@ -432,12 +456,9 @@ def sprepost(A, B): from .cy.qobjevo import QobjEvo if (isinstance(A, QobjEvo) or isinstance(B, QobjEvo)): return spre(A) * spost(B) - dims = [[_drop_projected_dims(A.dims[0]), - _drop_projected_dims(B.dims[1])], - [_drop_projected_dims(A.dims[1]), - _drop_projected_dims(B.dims[0])]] + return Qobj(_data.kron_transpose(B.data, A.data), - dims=dims, + dims=Dimensions.from_prepost(A._dims, B._dims), superrep='super', isherm=A._isherm and B._isherm, copy=False) diff --git a/qutip/core/transformation.py b/qutip/core/transformation.py new file mode 100644 index 0000000000..ef24674f9e --- /dev/null +++ b/qutip/core/transformation.py @@ -0,0 +1,70 @@ + +from .qobj import Qobj +from qutip.typing import LayerType +from . import data as _data +import qutip +from .dimensions import flatten +import numbers + +__all__ = ["ptrace"] + + +def ptrace(Q: Qobj, sel: int | list[int], dtype: LayerType = None) -> Qobj: + """ + Partial trace of the Qobj with selected components remaining. + + Parameters + ---------- + Q : :class:`.Qobj` + Composite quantum object. + sel : int/list + An ``int`` or ``list`` of components to keep after partial trace. + + Returns + ------- + oper : :class:`.Qobj` + Quantum object representing partial trace with selected components + remaining. + + dtype : type, str + The matrix format of output. + + Notes + ----- + This function is for legacy compatibility only. It is recommended to use + the ``ptrace()`` Qobj method. + """ + if not isinstance(Q, Qobj): + raise TypeError("Input is not a quantum object") + + try: + sel = sorted(sel) + except TypeError: + if not isinstance(sel, numbers.Integral): + raise TypeError( + "selection must be an integer or list of integers" + ) from None + sel = [sel] + if Q.isoperket: + dims = Q.dims[0] + data = qutip.vector_to_operator(Q).data + elif Q.isoperbra: + dims = Q.dims[1] + data = qutip.vector_to_operator(Q.dag()).data + elif Q.issuper or Q.isoper: + dims = Q.dims + data = Q.data + else: + dims = [Q.dims[0] if Q.isket else Q.dims[1]] * 2 + data = _data.project(Q.data) + if dims[0] != dims[1]: + raise ValueError("partial trace is not defined on non-square maps") + dims = flatten(dims[0]) + new_data = _data.ptrace(data, dims, sel, dtype=dtype) + new_dims = [[dims[x] for x in sel]] * 2 if sel else None + out = Qobj(new_data, dims=new_dims, copy=False) + if Q.isoperket: + return qutip.operator_to_vector(out) + if Q.isoperbra: + return qutip.operator_to_vector(out).dag() + return out diff --git a/qutip/solver/solver_base.py b/qutip/solver/solver_base.py index c41cce756a..75d4719871 100644 --- a/qutip/solver/solver_base.py +++ b/qutip/solver/solver_base.py @@ -101,7 +101,7 @@ def _prepare_state(self, state): } if state.isket: norm = state.norm() - elif state._dims.issquare: + elif state._dims.issquare and not state._dims.issuper: # Qobj.isoper does not differientiate between rectangular operators # and normal ones. norm = state.tr() diff --git a/qutip/tests/core/data/test_mathematics.py b/qutip/tests/core/data/test_mathematics.py index dc53bfbfeb..b64116b793 100644 --- a/qutip/tests/core/data/test_mathematics.py +++ b/qutip/tests/core/data/test_mathematics.py @@ -794,6 +794,24 @@ def op_numpy(self, left, right): ] +class TestMatmulDag(BinaryOpMixin): + def op_numpy(self, left, right): + return np.matmul(left, right) + + shapes = shapes_binary_matmul() + bad_shapes = shapes_binary_bad_matmul() + specialisations = [ + pytest.param( + lambda l, r: data.matmul_dag_data(l, r.adjoint()), + CSR, CSR, CSR + ), + pytest.param( + lambda l, r: data.matmul_dag_dense_csr_dense(l, r.adjoint()), + Dense, CSR, Dense + ), + ] + + class TestMultiply(BinaryOpMixin): def op_numpy(self, left, right): return left * right diff --git a/qutip/tests/core/test_gates.py b/qutip/tests/core/test_gates.py index 27f01cdf9d..55962f2188 100644 --- a/qutip/tests/core/test_gates.py +++ b/qutip/tests/core/test_gates.py @@ -4,6 +4,15 @@ from qutip.core import gates +def _calculate_isunitary(oper): + """ + Checks whether qobj is a unitary matrix + """ + if not oper.isoper or oper._data.shape[0] != oper._data.shape[1]: + return False + return oper @ oper.dag() == qutip.qeye_like(oper) + + def _infidelity(a, b): """Infidelity between two kets.""" return 1 - abs(a.overlap(b)) @@ -193,7 +202,7 @@ def test_metadata(gate_func, args, alias): dtype = qutip.data.to.parse(alias) assert isinstance(gate.data, dtype) assert gate._isherm == qutip.data.isherm(gate.data) - assert gate._isunitary == gate._calculate_isunitary() + assert gate._isunitary == _calculate_isunitary(gate) with qutip.CoreOptions(default_dtype=alias): gate = gate_func(*args) assert isinstance(gate.data, dtype) diff --git a/qutip/tests/core/test_metrics.py b/qutip/tests/core/test_metrics.py index def3e6b1a9..f8366484c3 100644 --- a/qutip/tests/core/test_metrics.py +++ b/qutip/tests/core/test_metrics.py @@ -205,21 +205,22 @@ def test_average_gate_fidelity_against_legacy_implementation(self): """ def agf_pre_50(oper, target=None): kraus = to_kraus(oper) - d = kraus[0].shape[0] + shape = kraus._pre_dims().shape - if kraus[0].shape[1] != d: + if shape[0] != shape[1]: return TypeError( "Average gate fidelity only implemented for square " "superoperators.") + d = shape[0] if target is None: return ( - (d + np.sum([np.abs(A_k.tr()) ** 2 for A_k in kraus])) / + (d + np.sum([np.abs(A_k.tr()) ** 2 for A_k in kraus._kmap])) / (d * d + d) ) return ( (d + np.sum([ - np.abs((A_k * target.dag()).tr()) ** 2 for A_k in kraus + np.abs((A_k * target.dag()).tr()) ** 2 for A_k in kraus._kmap ])) / (d * d + d) ) diff --git a/qutip/tests/core/test_operators.py b/qutip/tests/core/test_operators.py index 477d7da8a4..b13e819465 100644 --- a/qutip/tests/core/test_operators.py +++ b/qutip/tests/core/test_operators.py @@ -4,6 +4,16 @@ import qutip +def _calculate_isunitary(oper): + """ + Checks whether qobj is a unitary matrix + """ + if not oper.isoper or oper._data.shape[0] != oper._data.shape[1]: + return False + id_ = oper @ oper.dag() + return id_ == qutip.qeye_like(id_) + + N = 5 @@ -156,7 +166,7 @@ def test_position(): expected = (np.diag((np.arange(1, N) / 2)**0.5, k=-1) + np.diag((np.arange(1, N) / 2)**0.5, k=1)) np.testing.assert_allclose(operator.full(), expected) - assert operator._isherm == True + assert operator.isherm == True def test_momentum(): @@ -164,7 +174,7 @@ def test_momentum(): expected = (np.diag((np.arange(1, N) / 2)**0.5, k=-1) - np.diag((np.arange(1, N) / 2)**0.5, k=1)) * 1j np.testing.assert_allclose(operator.full(), expected) - assert operator._isherm == True + assert operator.isherm == True def test_squeeze(): @@ -252,7 +262,7 @@ def _check_meta(object, dtype): return assert isinstance(object.data, dtype) assert object._isherm == qutip.data.isherm(object.data) - assert object._isunitary == object._calculate_isunitary() + assert object._isunitary == _calculate_isunitary(object) # random object accept `str` and base.Data diff --git a/qutip/tests/core/test_qobj.py b/qutip/tests/core/test_qobj.py index 6e3ec30ee9..bce97affff 100644 --- a/qutip/tests/core/test_qobj.py +++ b/qutip/tests/core/test_qobj.py @@ -30,7 +30,7 @@ def assert_hermicity(oper, hermicity): # Check the cached isherm, if any exists. assert oper.isherm == hermicity # Force a reset of the cached value for isherm. - oper._isherm = None + oper.isherm = None # Force a recalculation of isherm. assert oper.isherm == hermicity @@ -162,7 +162,7 @@ def assert_unitarity(oper, unitarity): assert oper.isunitary == unitarity # Force a reset of the cached value for isunitary. - oper._isunitary = None + oper.isunitary = None # Force a recalculation of isunitary. assert oper.isunitary == unitarity @@ -261,7 +261,7 @@ def test_QobjAddition(): q4 = q1 + q2 q4_isherm = q4.isherm - q4._isherm = None # clear cached values + q4.isherm = None # clear cached values assert q4_isherm == q4.isherm # check elementwise addition/subtraction diff --git a/qutip/tests/core/test_superop_reps.py b/qutip/tests/core/test_superop_reps.py index 7190d0fe8c..661300651f 100644 --- a/qutip/tests/core/test_superop_reps.py +++ b/qutip/tests/core/test_superop_reps.py @@ -70,7 +70,8 @@ def test_SuperChoiSuper(self, superoperator): # Assert both that the result is close to expected, and has the right # type. - assert (test_supe - superoperator).norm() < tol + with CoreOptions(atol=tol): + assert test_supe == superoperator assert choi_matrix.type == "super" and choi_matrix.superrep == "choi" assert test_supe.type == "super" and test_supe.superrep == "super" @@ -90,7 +91,8 @@ def test_SuperChoiChiSuper(self, dimension): # Assert both that the result is close to expected, and has the right # type. - assert (test_supe - superoperator).norm() < tol + with CoreOptions(atol=tol): + assert test_supe == superoperator assert choi_matrix.type == "super" and choi_matrix.superrep == "choi" assert chi_matrix.type == "super" and chi_matrix.superrep == "chi" assert test_supe.type == "super" and test_supe.superrep == "super" @@ -105,7 +107,8 @@ def test_ChoiKrausChoi(self, superoperator): # Assert both that the result is close to expected, and has the right # type. - assert (test_choi - choi_matrix).norm() < tol + with CoreOptions(atol=tol): + assert test_choi == choi_matrix assert choi_matrix.type == "super" and choi_matrix.superrep == "choi" assert test_choi.type == "super" and test_choi.superrep == "choi" @@ -129,7 +132,8 @@ def test_NonSquareKrausSuperChoi(self): assert super.type == "super" and super.superrep == "super" assert_kraus_equivalence(op1[0], kraus, tol=1e-8) assert_kraus_equivalence(op2[0], kraus, tol=1e-8) - assert abs((op3 - super).norm()) < 1e-8 + with CoreOptions(atol=1e-8): + assert op3 == super def test_NeglectSmallKraus(self): """ @@ -233,8 +237,8 @@ def test_stinespring_cp(self, dimension): """ superop = rand_super_bcsz(dimension) A, B = to_stinespring(superop) - - assert (A - B).norm() < tol + with CoreOptions(atol=tol): + A == B @pytest.mark.repeat(3) def test_stinespring_agrees(self, dimension): @@ -255,7 +259,8 @@ def test_stinespring_agrees(self, dimension): # ptraced! q2 = (A * state * B.dag()).ptrace((0,)) - assert (q1 - q2).norm('tr') <= tol + with CoreOptions(atol=tol): + q1 == q2 def test_stinespring_dims(self, dimension): """ @@ -272,9 +277,9 @@ def test_chi_choi_roundtrip(self, dimension): superop = rand_super_bcsz(dimension) superop = to_chi(superop) rt_superop = to_chi(to_choi(superop)) - dif = (rt_superop - superop).norm() - assert dif == pytest.approx(0, abs=1e-7) + with CoreOptions(atol=1e-7): + assert rt_superop == superop assert rt_superop.type == superop.type assert rt_superop.dims == superop.dims @@ -312,4 +317,5 @@ def test_chi_known(self, superop, chi_expected): chiq = Qobj( chi_expected, dims=[[[2], [2]], [[2], [2]]], superrep='chi', ) - assert (chi_actual - chiq).norm() < tol + with CoreOptions(atol=tol): + assert chi_actual == chiq diff --git a/qutip/tests/solver/test_transfertensor.py b/qutip/tests/solver/test_transfertensor.py index d0f1061256..3154a097ed 100644 --- a/qutip/tests/solver/test_transfertensor.py +++ b/qutip/tests/solver/test_transfertensor.py @@ -27,8 +27,11 @@ def test_ttmsolve_jc_model(call): psi0c = qutip.basis(N, 0) rho0c = qutip.ket2dm(psi0c) rho0 = qutip.tensor(rho0a, rho0c) - superrho0cav = qutip.sprepost( - qutip.tensor(qutip.qeye(2), psi0c), qutip.tensor(qutip.qeye(2), psi0c.dag()) + superrho0cav = qutip.tensor_contract( + qutip.sprepost( + qutip.tensor(qutip.qeye(2), psi0c), + qutip.tensor(qutip.qeye(2), psi0c.dag()) + ), [5, 7] ) # calculate exact solution using mesolve