1# This file is part of QuTiP: Quantum Toolbox in Python.
2#
3#    Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
4#    All rights reserved.
5#
6#    Redistribution and use in source and binary forms, with or without
7#    modification, are permitted provided that the following conditions are
8#    met:
9#
10#    1. Redistributions of source code must retain the above copyright notice,
11#       this list of conditions and the following disclaimer.
12#
13#    2. Redistributions in binary form must reproduce the above copyright
14#       notice, this list of conditions and the following disclaimer in the
15#       documentation and/or other materials provided with the distribution.
16#
17#    3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
18#       of its contributors may be used to endorse or promote products derived
19#       from this software without specific prior written permission.
20#
21#    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22#    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23#    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
24#    PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25#    HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26#    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27#    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28#    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29#    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30#    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31#    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32###############################################################################
33"""The Quantum Object (Qobj) class, for representing quantum states and
34operators, and related functions.
35"""
36
37__all__ = ['Qobj', 'qobj_list_evaluate', 'ptrace', 'dag', 'isequal',
38           'issuper', 'isoper', 'isoperket', 'isoperbra', 'isket', 'isbra',
39           'isherm', 'shape', 'dims']
40
41import warnings
42import types
43import numbers
44
45try:
46    import builtins
47except ImportError:
48    import __builtin__ as builtins
49
50# import math functions from numpy.math: required for td string evaluation
51from numpy import (arccos, arccosh, arcsin, arcsinh, arctan, arctan2, arctanh,
52                   ceil, copysign, cos, cosh, degrees, e, exp, expm1, fabs,
53                   floor, fmod, frexp, hypot, isinf, isnan, ldexp, log, log10,
54                   log1p, modf, pi, radians, sin, sinh, sqrt, tan, tanh)
55
56import numpy as np
57import scipy.sparse as sp
58import qutip.settings as settings
59from qutip import __version__
60from qutip.fastsparse import fast_csr_matrix, fast_identity
61from qutip.cy.ptrace import _ptrace
62from qutip.permute import _permute
63from qutip.sparse import (
64    sp_eigs, sp_expm, sp_fro_norm, sp_max_norm, sp_one_norm, sp_L2_norm,
65)
66from qutip.dimensions import (
67    type_from_dims, enumerate_flat, collapse_dims_super,
68)
69from qutip.cy.spmath import (
70    zcsr_transpose, zcsr_adjoint, zcsr_isherm, zcsr_trace, zcsr_proj,
71    zcsr_inner,
72)
73from qutip.cy.spmatfuncs import zcsr_mat_elem
74from qutip.cy.sparse_utils import cy_tidyup
75import sys
76if sys.version_info.major >= 3:
77    from itertools import zip_longest
78elif sys.version_info.major < 3:
79    from itertools import izip_longest
80    zip_longest = izip_longest
81
82# OPENMP stuff
83from qutip.cy.openmp.utilities import use_openmp
84if settings.has_openmp:
85    from qutip.cy.openmp.omp_sparse_utils import omp_tidyup
86
87
88class Qobj(object):
89    """A class for representing quantum objects, such as quantum operators
90    and states.
91
92    The Qobj class is the QuTiP representation of quantum operators and state
93    vectors. This class also implements math operations +,-,* between Qobj
94    instances (and / by a C-number), as well as a collection of common
95    operator/state operations.  The Qobj constructor optionally takes a
96    dimension ``list`` and/or shape ``list`` as arguments.
97
98    Parameters
99    ----------
100    inpt : array_like
101        Data for vector/matrix representation of the quantum object.
102    dims : list
103        Dimensions of object used for tensor products.
104    shape : list
105        Shape of underlying data structure (matrix shape).
106    copy : bool
107        Flag specifying whether Qobj should get a copy of the
108        input data, or use the original.
109    fast : bool
110        Flag for fast qobj creation when running ode solvers.
111        This parameter is used internally only.
112
113
114    Attributes
115    ----------
116    data : array_like
117        Sparse matrix characterizing the quantum object.
118    dims : list
119        List of dimensions keeping track of the tensor structure.
120    shape : list
121        Shape of the underlying `data` array.
122    type : str
123        Type of quantum object: 'bra', 'ket', 'oper', 'operator-ket',
124        'operator-bra', or 'super'.
125    superrep : str
126        Representation used if `type` is 'super'. One of 'super'
127        (Liouville form) or 'choi' (Choi matrix with tr = dimension).
128    isherm : bool
129        Indicates if quantum object represents Hermitian operator.
130    isunitary : bool
131        Indictaes if quantum object represents unitary operator.
132    iscp : bool
133        Indicates if the quantum object represents a map, and if that map is
134        completely positive (CP).
135    ishp : bool
136        Indicates if the quantum object represents a map, and if that map is
137        hermicity preserving (HP).
138    istp : bool
139        Indicates if the quantum object represents a map, and if that map is
140        trace preserving (TP).
141    iscptp : bool
142        Indicates if the quantum object represents a map that is completely
143        positive and trace preserving (CPTP).
144    isket : bool
145        Indicates if the quantum object represents a ket.
146    isbra : bool
147        Indicates if the quantum object represents a bra.
148    isoper : bool
149        Indicates if the quantum object represents an operator.
150    issuper : bool
151        Indicates if the quantum object represents a superoperator.
152    isoperket : bool
153        Indicates if the quantum object represents an operator in column vector
154        form.
155    isoperbra : bool
156        Indicates if the quantum object represents an operator in row vector
157        form.
158
159    Methods
160    -------
161    copy()
162        Create copy of Qobj
163    conj()
164        Conjugate of quantum object.
165    cosm()
166        Cosine of quantum object.
167    dag()
168        Adjoint (dagger) of quantum object.
169    dnorm()
170        Diamond norm of quantum operator.
171    dual_chan()
172        Dual channel of quantum object representing a CP map.
173    eigenenergies(sparse=False, sort='low', eigvals=0, tol=0, maxiter=100000)
174        Returns eigenenergies (eigenvalues) of a quantum object.
175    eigenstates(sparse=False, sort='low', eigvals=0, tol=0, maxiter=100000)
176        Returns eigenenergies and eigenstates of quantum object.
177    expm()
178        Matrix exponential of quantum object.
179    full(order='C')
180        Returns dense array of quantum object `data` attribute.
181    groundstate(sparse=False, tol=0, maxiter=100000)
182        Returns eigenvalue and eigenket for the groundstate of a quantum
183        object.
184    inv()
185        Return a Qobj corresponding to the matrix inverse of the operator.
186    matrix_element(bra, ket)
187        Returns the matrix element of operator between `bra` and `ket` vectors.
188    norm(norm='tr', sparse=False, tol=0, maxiter=100000)
189        Returns norm of a ket or an operator.
190    permute(order)
191        Returns composite qobj with indices reordered.
192    proj()
193        Computes the projector for a ket or bra vector.
194    ptrace(sel)
195        Returns quantum object for selected dimensions after performing
196        partial trace.
197    sinm()
198        Sine of quantum object.
199    sqrtm()
200        Matrix square root of quantum object.
201    tidyup(atol=1e-12)
202        Removes small elements from quantum object.
203    tr()
204        Trace of quantum object.
205    trans()
206        Transpose of quantum object.
207    transform(inpt, inverse=False)
208        Performs a basis transformation defined by `inpt` matrix.
209    trunc_neg(method='clip')
210        Removes negative eigenvalues and returns a new Qobj that is
211        a valid density operator.
212    unit(norm='tr', sparse=False, tol=0, maxiter=100000)
213        Returns normalized quantum object.
214
215    """
216    __array_priority__ = 100  # sets Qobj priority above numpy arrays
217    # Disable ufuncs from acting directly on Qobj. This is necessary because we
218    # define __array__.
219    __array_ufunc__ = None
220
221    def __init__(self, inpt=None, dims=None, shape=None,
222                 type=None, isherm=None, copy=True,
223                 fast=False, superrep=None, isunitary=None):
224        """
225        Qobj constructor.
226        """
227        self._isherm = isherm
228        self._type = type
229        self.superrep = superrep
230        self._isunitary = isunitary
231
232        if fast == 'mc':
233            # fast Qobj construction for use in mcsolve with ket output
234            self._data = inpt
235            self.dims = dims
236            self._isherm = False
237            return
238
239        if fast == 'mc-dm':
240            # fast Qobj construction for use in mcsolve with dm output
241            self._data = inpt
242            self.dims = dims
243            self._isherm = True
244            return
245
246        if isinstance(inpt, Qobj):
247            # if input is already Qobj then return identical copy
248
249            self._data = fast_csr_matrix(
250                (inpt.data.data, inpt.data.indices, inpt.data.indptr),
251                shape=inpt.shape, copy=copy,
252            )
253
254            if dims is None:
255                # Dimensions of quantum object used for keeping track of tensor
256                # components
257                self.dims = inpt.dims
258            else:
259                self.dims = dims
260
261            self.superrep = inpt.superrep
262            self._isunitary = inpt._isunitary
263
264        elif inpt is None:
265            # initialize an empty Qobj with correct dimensions and shape
266
267            if dims is not None:
268                N, M = np.prod(dims[0]), np.prod(dims[1])
269                self.dims = dims
270
271            elif shape is not None:
272                N, M = shape
273                self.dims = [[N], [M]]
274
275            else:
276                N, M = 1, 1
277                self.dims = [[N], [M]]
278
279            self._data = fast_csr_matrix(shape=(N, M))
280
281        elif isinstance(inpt, list) or isinstance(inpt, tuple):
282            # case where input is a list
283            data = np.array(inpt)
284            if len(data.shape) == 1:
285                # if list has only one dimension (i.e [5,4])
286                data = data.transpose()
287
288            _tmp = sp.csr_matrix(data, dtype=complex)
289            self._data = fast_csr_matrix(
290                (_tmp.data, _tmp.indices, _tmp.indptr),
291                shape=_tmp.shape,
292            )
293            if dims is None:
294                self.dims = [[int(data.shape[0])], [int(data.shape[1])]]
295            else:
296                self.dims = dims
297
298        elif isinstance(inpt, np.ndarray) or sp.issparse(inpt):
299            # case where input is array or sparse
300            if inpt.ndim == 1:
301                inpt = inpt[:, np.newaxis]
302
303            do_copy = copy
304            if not isinstance(inpt, fast_csr_matrix):
305                _tmp = sp.csr_matrix(inpt, dtype=complex, copy=do_copy)
306                _tmp.sort_indices()  # Make sure indices are sorted.
307                do_copy = 0
308            else:
309                _tmp = inpt
310            self._data = fast_csr_matrix(
311                (_tmp.data, _tmp.indices, _tmp.indptr),
312                shape=_tmp.shape,
313                copy=do_copy,
314            )
315
316            if dims is None:
317                self.dims = [[int(inpt.shape[0])], [int(inpt.shape[1])]]
318            else:
319                self.dims = dims
320
321        elif isinstance(inpt, (int, float, complex,
322                               np.integer, np.floating, np.complexfloating)):
323            # if input is int, float, or complex then convert to array
324            _tmp = sp.csr_matrix([[inpt]], dtype=complex)
325            self._data = fast_csr_matrix(
326                (_tmp.data, _tmp.indices, _tmp.indptr),
327                shape=_tmp.shape,
328            )
329            if dims is None:
330                self.dims = [[1], [1]]
331            else:
332                self.dims = dims
333
334        else:
335            warnings.warn("Initializing Qobj from unsupported type: %s" %
336                          builtins.type(inpt))
337            inpt = np.array([[0]])
338            _tmp = sp.csr_matrix(inpt, dtype=complex, copy=copy)
339            self._data = fast_csr_matrix(
340                (_tmp.data, _tmp.indices, _tmp.indptr),
341                shape=_tmp.shape,
342            )
343            self.dims = [[int(inpt.shape[0])], [int(inpt.shape[1])]]
344
345        if type == 'super':
346            # Type is not super, i.e. dims not explicitly passed, but oper-like
347            # shape.
348            if dims is None and self.shape[0] == self.shape[1]:
349                sub_shape = np.sqrt(self.shape[0])
350                # check if root of shape is int
351                if (sub_shape % 1) != 0:
352                    raise Exception('Invalid shape for a super operator.')
353                else:
354                    sub_shape = int(sub_shape)
355                    self.dims = [[[sub_shape], [sub_shape]]]*2
356
357        if superrep:
358            self.superrep = superrep
359        else:
360            if self.type == 'super' and self.superrep is None:
361                self.superrep = 'super'
362
363        # clear type cache
364        self._type = None
365
366    def copy(self):
367        """Create identical copy"""
368        return Qobj(inpt=self)
369
370    def get_data(self):
371        return self._data
372
373    # Here we perfrom a check of the csr matrix type during setting of Q.data
374    def set_data(self, data):
375        if not isinstance(data, fast_csr_matrix):
376            raise TypeError('Qobj data must be in fast_csr format.')
377        else:
378            self._data = data
379    data = property(get_data, set_data)
380
381    def __add__(self, other):
382        """
383        ADDITION with Qobj on LEFT [ ex. Qobj+4 ]
384        """
385        self._isunitary = None
386
387        if isinstance(other, eseries):
388            return other.__radd__(self)
389
390        if not isinstance(other, Qobj):
391            if isinstance(other, (int, float, complex, np.integer, np.floating,
392                          np.complexfloating, np.ndarray, list, tuple)) \
393                          or sp.issparse(other):
394                other = Qobj(other)
395            else:
396                return NotImplemented
397
398        if np.prod(other.shape) == 1 and np.prod(self.shape) != 1:
399            # case for scalar quantum object
400            dat = other.data[0, 0]
401            if dat == 0:
402                return self
403
404            out = Qobj()
405
406            if self.type in ['oper', 'super']:
407                out.data = self.data + dat * fast_identity(
408                    self.shape[0])
409            else:
410                out.data = self.data
411                out.data.data = out.data.data + dat
412
413            out.dims = self.dims
414
415            if settings.auto_tidyup:
416                out.tidyup()
417
418            if isinstance(dat, (int, float)):
419                out._isherm = self._isherm
420            else:
421                # We use _isherm here to prevent recalculating on self and
422                # other, relying on that bool(None) == False.
423                out._isherm = (self._isherm and other._isherm) or out.isherm
424
425            out.superrep = self.superrep
426
427            return out
428
429        elif np.prod(self.shape) == 1 and np.prod(other.shape) != 1:
430            # case for scalar quantum object
431            dat = self.data[0, 0]
432            if dat == 0:
433                return other
434
435            out = Qobj()
436            if other.type in ['oper', 'super']:
437                out.data = dat * fast_identity(other.shape[0]) + other.data
438            else:
439                out.data = other.data
440                out.data.data = out.data.data + dat
441            out.dims = other.dims
442
443            if settings.auto_tidyup:
444                out.tidyup()
445
446            if isinstance(dat, complex):
447                out._isherm = out.isherm
448            else:
449                out._isherm = self._isherm
450
451            out.superrep = self.superrep
452
453            return out
454
455        elif self.dims != other.dims:
456            raise TypeError('Incompatible quantum object dimensions')
457
458        elif self.shape != other.shape:
459            raise TypeError('Matrix shapes do not match')
460
461        else:  # case for matching quantum objects
462            out = Qobj()
463            out.data = self.data + other.data
464            out.dims = self.dims
465            if settings.auto_tidyup:
466                out.tidyup()
467
468            if self.type in ['ket', 'bra', 'operator-ket', 'operator-bra']:
469                out._isherm = False
470            elif self._isherm is None or other._isherm is None:
471                out._isherm = out.isherm
472            elif not self._isherm and not other._isherm:
473                out._isherm = out.isherm
474            else:
475                out._isherm = self._isherm and other._isherm
476
477            if self.superrep and other.superrep:
478                if self.superrep != other.superrep:
479                    msg = ("Adding superoperators with different " +
480                           "representations")
481                    warnings.warn(msg)
482
483                out.superrep = self.superrep
484
485            return out
486
487    def __radd__(self, other):
488        """
489        ADDITION with Qobj on RIGHT [ ex. 4+Qobj ]
490        """
491        return self + other
492
493    def __sub__(self, other):
494        """
495        SUBTRACTION with Qobj on LEFT [ ex. Qobj-4 ]
496        """
497        return self + (-other)
498
499    def __rsub__(self, other):
500        """
501        SUBTRACTION with Qobj on RIGHT [ ex. 4-Qobj ]
502        """
503        return (-self) + other
504
505    def __mul__(self, other):
506        """
507        MULTIPLICATION with Qobj on LEFT [ ex. Qobj*4 ]
508        """
509        self._isunitary = None
510
511        if isinstance(other, Qobj):
512            if self.dims[1] == other.dims[0]:
513                out = Qobj()
514                out.data = self.data * other.data
515                dims = [self.dims[0], other.dims[1]]
516                out.dims = dims
517                if settings.auto_tidyup:
518                    out.tidyup()
519                if (settings.auto_tidyup_dims
520                        and not isinstance(dims[0][0], list)
521                        and not isinstance(dims[1][0], list)):
522                    # If neither left or right is a superoperator,
523                    # we should implicitly partial trace over
524                    # matching dimensions of 1.
525                    # Using izip_longest allows for the left and right dims
526                    # to have uneven length (non-square Qobjs).
527                    # We use None as padding so that it doesn't match anything,
528                    # and will never cause a partial trace on the other side.
529                    mask = [l == r == 1
530                            for l, r in zip_longest(dims[0], dims[1],
531                                                    fillvalue=None)]
532                    # To ensure that there are still any dimensions left, we
533                    # use max() to add a dimensions list of [1] if all matching
534                    # dims are traced out of that side.
535                    out.dims = [max([1],
536                                    [dim for dim, m in zip(dims[0], mask)
537                                    if not m]),
538                                max([1],
539                                    [dim for dim, m in zip(dims[1], mask)
540                                    if not m])]
541
542                else:
543                    out.dims = dims
544
545                out._isherm = None
546
547                if self.superrep and other.superrep:
548                    if self.superrep != other.superrep:
549                        msg = ("Multiplying superoperators with different " +
550                               "representations")
551                        warnings.warn(msg)
552
553                    out.superrep = self.superrep
554
555                return out
556
557            elif np.prod(self.shape) == 1:
558                out = Qobj(other)
559                out.data *= self.data[0, 0]
560                out.superrep = other.superrep
561                return out.tidyup() if settings.auto_tidyup else out
562
563            elif np.prod(other.shape) == 1:
564                out = Qobj(self)
565                out.data *= other.data[0, 0]
566                out.superrep = self.superrep
567                return out.tidyup() if settings.auto_tidyup else out
568
569            elif self.shape[1] == 1 and other.shape[0] == 1:
570                out = Qobj()
571                out.data = self.data * other.data
572                out.dims = [self.dims[0], other.dims[1]]
573                return out.tidyup() if settings.auto_tidyup else out
574
575            else:
576                raise TypeError("Incompatible Qobj shapes")
577
578        elif isinstance(other, np.ndarray):
579            if other.dtype == 'object':
580                out = np.empty(other.shape, dtype=object)
581                for i, item in enumerate(other.flat):
582                    out.flat[i] = self * item
583                return out
584            else:
585                return self.data * other
586
587        elif isinstance(other, list):
588            # if other is a list, do element-wise multiplication
589            out = np.empty(len(other), dtype=object)
590            out[:] = [self * item for item in other]
591            return out
592
593        elif isinstance(other, eseries):
594            return other.__rmul__(self)
595
596        elif isinstance(other, numbers.Number):
597            out = Qobj()
598            out.data = self.data * other
599            out.dims = self.dims
600            out.superrep = self.superrep
601            if settings.auto_tidyup:
602                out.tidyup()
603            if isinstance(other, complex):
604                out._isherm = out.isherm
605            else:
606                out._isherm = self._isherm
607
608            return out
609
610        else:
611            return NotImplemented
612
613    def __rmul__(self, other):
614        """
615        MULTIPLICATION with Qobj on RIGHT [ ex. 4*Qobj ]
616        """
617        if isinstance(other, np.ndarray):
618            if other.dtype == 'object':
619                out = np.empty(other.shape, dtype=object)
620                for i, item in enumerate(other.flat):
621                    out.flat[i] = item * self
622                return out
623            else:
624                return other * self.data
625
626        elif isinstance(other, list):
627            # if other is a list, do element-wise multiplication
628            out = np.empty(len(other), dtype=object)
629            out[:] = [item * self for item in other]
630            return out
631
632        elif isinstance(other, eseries):
633            return other.__mul__(self)
634
635        elif isinstance(other, numbers.Number):
636            out = Qobj()
637            out.data = other * self.data
638            out.dims = self.dims
639            out.superrep = self.superrep
640            if settings.auto_tidyup:
641                out.tidyup()
642            if isinstance(other, complex):
643                out._isherm = out.isherm
644            else:
645                out._isherm = self._isherm
646
647            return out
648
649        else:
650            raise TypeError("Incompatible object for multiplication")
651
652    def __truediv__(self, other):
653        return self.__div__(other)
654
655    def __div__(self, other):
656        """
657        DIVISION (by numbers only)
658        """
659        if isinstance(other, Qobj):  # if both are quantum objects
660            raise TypeError("Incompatible Qobj shapes " +
661                            "[division with Qobj not implemented]")
662
663        if isinstance(other, (int, float, complex,
664                              np.integer, np.floating, np.complexfloating)):
665            out = Qobj()
666            out.data = self.data / other
667            out.dims = self.dims
668            if settings.auto_tidyup:
669                out.tidyup()
670            if isinstance(other, complex):
671                out._isherm = out.isherm
672            else:
673                out._isherm = self._isherm
674
675            out.superrep = self.superrep
676
677            return out
678
679        else:
680            raise TypeError("Incompatible object for division")
681
682    def __neg__(self):
683        """
684        NEGATION operation.
685        """
686        out = Qobj()
687        out.data = -self.data
688        out.dims = self.dims
689        out.superrep = self.superrep
690        if settings.auto_tidyup:
691            out.tidyup()
692        out._isherm = self._isherm
693        out._isunitary = self._isunitary
694        return out
695
696    def __getitem__(self, ind):
697        """
698        GET qobj elements.
699        """
700        out = self.data[ind]
701        if sp.issparse(out):
702            return out.toarray()
703        else:
704            return out
705
706    def __eq__(self, other):
707        """
708        EQUALITY operator.
709        """
710        if (isinstance(other, Qobj) and
711                self.dims == other.dims and
712                not np.any(np.abs((self.data - other.data).data) >
713                           settings.atol)):
714            return True
715        else:
716            return False
717
718    def __ne__(self, other):
719        """
720        INEQUALITY operator.
721        """
722        return not (self == other)
723
724    def __pow__(self, n, m=None):  # calculates powers of Qobj
725        """
726        POWER operation.
727        """
728        if self.type not in ['oper', 'super']:
729            raise Exception("Raising a qobj to some power works only for " +
730                            "operators and super-operators (square matrices).")
731
732        if m is not None:
733            raise NotImplementedError("modulo is not implemented for Qobj")
734
735        try:
736            data = self.data ** n
737        except (TypeError, ValueError):
738            raise ValueError('Invalid choice of exponent.')
739        out = Qobj(data, dims=self.dims)
740        out.superrep = self.superrep
741        return out.tidyup() if settings.auto_tidyup else out
742
743    def __abs__(self):
744        return abs(self.data)
745
746    def __str__(self):
747        s = ""
748        t = self.type
749        shape = self.shape
750        if self.type in ['oper', 'super']:
751            s += ("Quantum object: " +
752                  "dims = " + str(self.dims) +
753                  ", shape = " + str(shape) +
754                  ", type = " + t +
755                  ", isherm = " + str(self.isherm) +
756                  (
757                      ", superrep = {0.superrep}".format(self)
758                      if t == "super" and self.superrep != "super"
759                      else ""
760                  ) + "\n")
761        else:
762            s += ("Quantum object: " +
763                  "dims = " + str(self.dims) +
764                  ", shape = " + str(shape) +
765                  ", type = " + t + "\n")
766        s += "Qobj data =\n"
767
768        if shape[0] > 10000 or shape[1] > 10000:
769            # if the system is huge, don't attempt to convert to a
770            # dense matrix and then to string, because it is pointless
771            # and is likely going to produce memory errors. Instead print the
772            # sparse data string representation
773            s += str(self.data)
774
775        elif all(np.imag(self.data.data) == 0):
776            s += str(np.real(self.full()))
777
778        else:
779            s += str(self.full())
780
781        return s
782
783    def __repr__(self):
784        # give complete information on Qobj without print statement in
785        # command-line we cant realistically serialize a Qobj into a string,
786        # so we simply return the informal __str__ representation instead.)
787        return self.__str__()
788
789    def __call__(self, other):
790        """
791        Acts this Qobj on another Qobj either by left-multiplication,
792        or by vectorization and devectorization, as
793        appropriate.
794        """
795        if not isinstance(other, Qobj):
796            raise TypeError("Only defined for quantum objects.")
797
798        if self.type == "super":
799            if other.type == "ket":
800                other = qutip.states.ket2dm(other)
801
802            if other.type == "oper":
803                return qutip.superoperator.vector_to_operator(
804                    self * qutip.superoperator.operator_to_vector(other)
805                )
806            else:
807                raise TypeError("Can only act super on oper or ket.")
808
809        elif self.type == "oper":
810            if other.type == "ket":
811                return self * other
812            else:
813                raise TypeError("Can only act oper on ket.")
814
815    def __getstate__(self):
816        # defines what happens when Qobj object gets pickled
817        self.__dict__.update({'qutip_version': __version__[:5]})
818        return self.__dict__
819
820    def __setstate__(self, state):
821        # defines what happens when loading a pickled Qobj
822        if 'qutip_version' in state.keys():
823            del state['qutip_version']
824        (self.__dict__).update(state)
825
826    def _repr_latex_(self):
827        """
828        Generate a LaTeX representation of the Qobj instance. Can be used for
829        formatted output in ipython notebook.
830        """
831        t = self.type
832        shape = self.shape
833        s = r''
834        if self.type in ['oper', 'super']:
835            s += ("Quantum object: " +
836                  "dims = " + str(self.dims) +
837                  ", shape = " + str(shape) +
838                  ", type = " + t +
839                  ", isherm = " + str(self.isherm) +
840                  (
841                      ", superrep = {0.superrep}".format(self)
842                      if t == "super" and self.superrep != "super"
843                      else ""
844                  ))
845        else:
846            s += ("Quantum object: " +
847                  "dims = " + str(self.dims) +
848                  ", shape = " + str(shape) +
849                  ", type = " + t)
850
851        M, N = self.data.shape
852
853        s += r'\begin{equation*}\left(\begin{array}{*{11}c}'
854
855        def _format_float(value):
856            if value == 0.0:
857                return "0.0"
858            elif abs(value) > 1000.0 or abs(value) < 0.001:
859                return ("%.3e" % value).replace("e", r"\times10^{") + "}"
860            elif abs(value - int(value)) < 0.001:
861                return "%.1f" % value
862            else:
863                return "%.3f" % value
864
865        def _format_element(m, n, d):
866            s = " & " if n > 0 else ""
867            if type(d) == str:
868                return s + d
869            else:
870                if abs(np.imag(d)) < settings.atol:
871                    return s + _format_float(np.real(d))
872                elif abs(np.real(d)) < settings.atol:
873                    return s + _format_float(np.imag(d)) + "j"
874                else:
875                    s_re = _format_float(np.real(d))
876                    s_im = _format_float(np.imag(d))
877                    if np.imag(d) > 0.0:
878                        return (s + "(" + s_re + "+" + s_im + "j)")
879                    else:
880                        return (s + "(" + s_re + s_im + "j)")
881
882        if M > 10 and N > 10:
883            # truncated matrix output
884            for m in range(5):
885                for n in range(5):
886                    s += _format_element(m, n, self.data[m, n])
887                s += r' & \cdots'
888                for n in range(N - 5, N):
889                    s += _format_element(m, n, self.data[m, n])
890                s += r'\\'
891
892            for n in range(5):
893                s += _format_element(m, n, r'\vdots')
894            s += r' & \ddots'
895            for n in range(N - 5, N):
896                s += _format_element(m, n, r'\vdots')
897            s += r'\\'
898
899            for m in range(M - 5, M):
900                for n in range(5):
901                    s += _format_element(m, n, self.data[m, n])
902                s += r' & \cdots'
903                for n in range(N - 5, N):
904                    s += _format_element(m, n, self.data[m, n])
905                s += r'\\'
906
907        elif M > 10 and N <= 10:
908            # truncated vertically elongated matrix output
909            for m in range(5):
910                for n in range(N):
911                    s += _format_element(m, n, self.data[m, n])
912                s += r'\\'
913
914            for n in range(N):
915                s += _format_element(m, n, r'\vdots')
916            s += r'\\'
917
918            for m in range(M - 5, M):
919                for n in range(N):
920                    s += _format_element(m, n, self.data[m, n])
921                s += r'\\'
922
923        elif M <= 10 and N > 10:
924            # truncated horizontally elongated matrix output
925            for m in range(M):
926                for n in range(5):
927                    s += _format_element(m, n, self.data[m, n])
928                s += r' & \cdots'
929                for n in range(N - 5, N):
930                    s += _format_element(m, n, self.data[m, n])
931                s += r'\\'
932
933        else:
934            # full output
935            for m in range(M):
936                for n in range(N):
937                    s += _format_element(m, n, self.data[m, n])
938                s += r'\\'
939
940        s += r'\end{array}\right)\end{equation*}'
941        return s
942
943    def dag(self):
944        """Adjoint operator of quantum object.
945        """
946        out = Qobj()
947        out.data = zcsr_adjoint(self.data)
948        out.dims = [self.dims[1], self.dims[0]]
949        out._isherm = self._isherm
950        out.superrep = self.superrep
951        return out
952
953    def dual_chan(self):
954        """Dual channel of quantum object representing a completely positive
955        map.
956        """
957        # Uses the technique of Johnston and Kribs (arXiv:1102.0948), which
958        # is only valid for completely positive maps.
959        if not self.iscp:
960            raise ValueError("Dual channels are only implemented for CP maps.")
961        J = sr.to_choi(self)
962        tensor_idxs = enumerate_flat(J.dims)
963        J_dual = tensor.tensor_swap(J, *(
964                list(zip(tensor_idxs[0][1], tensor_idxs[0][0])) +
965                list(zip(tensor_idxs[1][1], tensor_idxs[1][0]))
966        )).trans()
967        J_dual.superrep = 'choi'
968        return J_dual
969
970    def conj(self):
971        """Conjugate operator of quantum object.
972        """
973        out = Qobj()
974        out.data = self.data.conj()
975        out.dims = [self.dims[0], self.dims[1]]
976        return out
977
978    def norm(self, norm=None, sparse=False, tol=0, maxiter=100000):
979        """Norm of a quantum object.
980
981        Default norm is L2-norm for kets and trace-norm for operators.
982        Other ket and operator norms may be specified using the `norm` and
983        argument.
984
985        Parameters
986        ----------
987        norm : str
988            Which norm to use for ket/bra vectors: L2 'l2', max norm 'max',
989            or for operators: trace 'tr', Frobius 'fro', one 'one', or max
990            'max'.
991
992        sparse : bool
993            Use sparse eigenvalue solver for trace norm.  Other norms are not
994            affected by this parameter.
995
996        tol : float
997            Tolerance for sparse solver (if used) for trace norm. The sparse
998            solver may not converge if the tolerance is set too low.
999
1000        maxiter : int
1001            Maximum number of iterations performed by sparse solver (if used)
1002            for trace norm.
1003
1004        Returns
1005        -------
1006        norm : float
1007            The requested norm of the operator or state quantum object.
1008
1009
1010        Notes
1011        -----
1012        The sparse eigensolver is much slower than the dense version.
1013        Use sparse only if memory requirements demand it.
1014
1015        """
1016        if self.type in ['oper', 'super']:
1017            if norm is None or norm == 'tr':
1018                _op = self.data * zcsr_adjoint(self.data)
1019                vals = sp_eigs(_op, True, vecs=False,
1020                               sparse=sparse, tol=tol, maxiter=maxiter)
1021                return np.sum(np.sqrt(np.abs(vals)))
1022            elif norm == 'fro':
1023                return sp_fro_norm(self.data)
1024            elif norm == 'one':
1025                return sp_one_norm(self.data)
1026            elif norm == 'max':
1027                return sp_max_norm(self.data)
1028            else:
1029                raise ValueError(
1030                    "For matrices, norm must be 'tr', 'fro', 'one', or 'max'.")
1031        else:
1032            if norm is None or norm == 'l2':
1033                return sp_L2_norm(self.data)
1034            elif norm == 'max':
1035                return sp_max_norm(self.data)
1036            else:
1037                raise ValueError("For vectors, norm must be 'l2', or 'max'.")
1038
1039    def proj(self):
1040        """Form the projector from a given ket or bra vector.
1041
1042        Parameters
1043        ----------
1044        Q : :class:`qutip.Qobj`
1045            Input bra or ket vector
1046
1047        Returns
1048        -------
1049        P : :class:`qutip.Qobj`
1050            Projection operator.
1051        """
1052        if self.isket:
1053            _out = zcsr_proj(self.data, 1)
1054            _dims = [self.dims[0], self.dims[0]]
1055        elif self.isbra:
1056            _out = zcsr_proj(self.data, 0)
1057            _dims = [self.dims[1], self.dims[1]]
1058        else:
1059            raise TypeError('Projector can only be formed from a bra or ket.')
1060
1061        return Qobj(_out, dims=_dims)
1062
1063    def tr(self):
1064        """Trace of a quantum object.
1065
1066        Returns
1067        -------
1068        trace : float
1069            Returns the trace of the quantum object.
1070
1071        """
1072        return zcsr_trace(self.data, self.isherm)
1073
1074    def purity(self):
1075        """Calculate purity of a quantum object.
1076
1077        Returns
1078        -------
1079        state_purity : float
1080            Returns the purity of a quantum object.
1081            For a pure state, the purity is 1.
1082            For a mixed state of dimension `d`, 1/d<=purity<1.
1083
1084        """
1085        rho = self
1086        if (rho.type == "super"):
1087            raise TypeError('Purity is defined on a density matrix or state.')
1088
1089        if rho.type == "ket" or rho.type == "bra":
1090            state_purity = (rho*rho.dag()).tr()
1091
1092        if rho.type == "oper":
1093            state_purity = (rho*rho).tr()
1094
1095        return state_purity
1096
1097    def full(self, order='C', squeeze=False):
1098        """Dense array from quantum object.
1099
1100        Parameters
1101        ----------
1102        order : str {'C', 'F'}
1103            Return array in C (default) or Fortran ordering.
1104        squeeze : bool {False, True}
1105            Squeeze output array.
1106
1107        Returns
1108        -------
1109        data : array
1110            Array of complex data from quantum objects `data` attribute.
1111        """
1112        if squeeze:
1113            return self.data.toarray(order=order).squeeze()
1114        else:
1115            return self.data.toarray(order=order)
1116
1117    def __array__(self, *arg, **kwarg):
1118        """Numpy array from Qobj
1119        For compatibility with np.array
1120        """
1121        return self.full()
1122
1123    def diag(self):
1124        """Diagonal elements of quantum object.
1125
1126        Returns
1127        -------
1128        diags : array
1129            Returns array of ``real`` values if operators is Hermitian,
1130            otherwise ``complex`` values are returned.
1131
1132        """
1133        out = self.data.diagonal()
1134        if np.any(np.imag(out) > settings.atol) or not self.isherm:
1135            return out
1136        else:
1137            return np.real(out)
1138
1139    def expm(self, method='dense'):
1140        """Matrix exponential of quantum operator.
1141
1142        Input operator must be square.
1143
1144        Parameters
1145        ----------
1146        method : str {'dense', 'sparse'}
1147            Use set method to use to calculate the matrix exponentiation. The
1148            available choices includes 'dense' and 'sparse'.  Since the
1149            exponential of a matrix is nearly always dense, method='dense'
1150            is set as default.s
1151
1152        Returns
1153        -------
1154        oper : :class:`qutip.Qobj`
1155            Exponentiated quantum operator.
1156
1157        Raises
1158        ------
1159        TypeError
1160            Quantum operator is not square.
1161
1162        """
1163        if self.dims[0][0] != self.dims[1][0]:
1164            raise TypeError('Invalid operand for matrix exponential')
1165
1166        if method == 'dense':
1167            F = sp_expm(self.data, sparse=False)
1168
1169        elif method == 'sparse':
1170            F = sp_expm(self.data, sparse=True)
1171
1172        else:
1173            raise ValueError("method must be 'dense' or 'sparse'.")
1174
1175        out = Qobj(F, dims=self.dims)
1176        return out.tidyup() if settings.auto_tidyup else out
1177
1178    def check_herm(self):
1179        """Check if the quantum object is hermitian.
1180
1181        Returns
1182        -------
1183        isherm : bool
1184            Returns the new value of isherm property.
1185        """
1186        self._isherm = None
1187        return self.isherm
1188
1189    def sqrtm(self, sparse=False, tol=0, maxiter=100000):
1190        """Sqrt of a quantum operator.
1191
1192        Operator must be square.
1193
1194        Parameters
1195        ----------
1196        sparse : bool
1197            Use sparse eigenvalue/vector solver.
1198        tol : float
1199            Tolerance used by sparse solver (0 = machine precision).
1200        maxiter : int
1201            Maximum number of iterations used by sparse solver.
1202
1203        Returns
1204        -------
1205        oper : :class:`qutip.Qobj`
1206            Matrix square root of operator.
1207
1208        Raises
1209        ------
1210        TypeError
1211            Quantum object is not square.
1212
1213        Notes
1214        -----
1215        The sparse eigensolver is much slower than the dense version.
1216        Use sparse only if memory requirements demand it.
1217
1218        """
1219        if self.dims[0][0] == self.dims[1][0]:
1220            evals, evecs = sp_eigs(self.data, self.isherm, sparse=sparse,
1221                                   tol=tol, maxiter=maxiter)
1222            numevals = len(evals)
1223            dV = sp.spdiags(np.sqrt(evals, dtype=complex), 0, numevals,
1224                            numevals, format='csr')
1225            if self.isherm:
1226                spDv = dV.dot(evecs.T.conj().T)
1227            else:
1228                spDv = dV.dot(np.linalg.inv(evecs.T))
1229
1230            out = Qobj(evecs.T.dot(spDv), dims=self.dims)
1231            return out.tidyup() if settings.auto_tidyup else out
1232
1233        else:
1234            raise TypeError('Invalid operand for matrix square root')
1235
1236    def cosm(self):
1237        """Cosine of a quantum operator.
1238
1239        Operator must be square.
1240
1241        Returns
1242        -------
1243        oper : :class:`qutip.Qobj`
1244            Matrix cosine of operator.
1245
1246        Raises
1247        ------
1248        TypeError
1249            Quantum object is not square.
1250
1251        Notes
1252        -----
1253        Uses the Q.expm() method.
1254
1255        """
1256        if self.dims[0][0] == self.dims[1][0]:
1257            return 0.5 * ((1j * self).expm() + (-1j * self).expm())
1258        else:
1259            raise TypeError('Invalid operand for matrix square root')
1260
1261    def sinm(self):
1262        """Sine of a quantum operator.
1263
1264        Operator must be square.
1265
1266        Returns
1267        -------
1268        oper : :class:`qutip.Qobj`
1269            Matrix sine of operator.
1270
1271        Raises
1272        ------
1273        TypeError
1274            Quantum object is not square.
1275
1276        Notes
1277        -----
1278        Uses the Q.expm() method.
1279
1280        """
1281        if self.dims[0][0] == self.dims[1][0]:
1282            return -0.5j * ((1j * self).expm() - (-1j * self).expm())
1283        else:
1284            raise TypeError('Invalid operand for matrix square root')
1285
1286    def inv(self, sparse=False):
1287        """Matrix inverse of a quantum operator
1288
1289        Operator must be square.
1290
1291        Returns
1292        -------
1293        oper : :class:`qutip.Qobj`
1294            Matrix inverse of operator.
1295
1296        Raises
1297        ------
1298        TypeError
1299            Quantum object is not square.
1300        """
1301        if self.shape[0] != self.shape[1]:
1302            raise TypeError('Invalid operand for matrix inverse')
1303        if sparse:
1304            inv_mat = sp.linalg.inv(sp.csc_matrix(self.data))
1305        else:
1306            inv_mat = np.linalg.inv(self.full())
1307        return Qobj(inv_mat, dims=self.dims[::-1])
1308
1309    def unit(self, inplace=False,
1310             norm=None, sparse=False,
1311             tol=0, maxiter=100000):
1312        """Operator or state normalized to unity.
1313
1314        Uses norm from Qobj.norm().
1315
1316        Parameters
1317        ----------
1318        inplace : bool
1319            Do an in-place normalization
1320        norm : str
1321            Requested norm for states / operators.
1322        sparse : bool
1323            Use sparse eigensolver for trace norm. Does not affect other norms.
1324        tol : float
1325            Tolerance used by sparse eigensolver.
1326        maxiter : int
1327            Number of maximum iterations performed by sparse eigensolver.
1328
1329        Returns
1330        -------
1331        oper : :class:`qutip.Qobj`
1332            Normalized quantum object if not in-place,
1333            else None.
1334
1335        """
1336        if inplace:
1337            nrm = self.norm(norm=norm, sparse=sparse,
1338                            tol=tol, maxiter=maxiter)
1339
1340            self.data /= nrm
1341        elif not inplace:
1342            out = self / self.norm(norm=norm, sparse=sparse,
1343                                   tol=tol, maxiter=maxiter)
1344            if settings.auto_tidyup:
1345                return out.tidyup()
1346            else:
1347                return out
1348        else:
1349            raise Exception('inplace kwarg must be bool.')
1350
1351    def ptrace(self, sel, sparse=None):
1352        """Partial trace of the quantum object.
1353
1354        Parameters
1355        ----------
1356        sel : int/list
1357            An ``int`` or ``list`` of components to keep after partial trace.
1358            The order is unimportant; no transposition will be done and the
1359            spaces will remain in the same order in the output.
1360
1361        Returns
1362        -------
1363        oper : :class:`qutip.Qobj`
1364            Quantum object representing partial trace with selected components
1365            remaining.
1366
1367        Notes
1368        -----
1369        This function is identical to the :func:`qutip.qobj.ptrace` function
1370        that has been deprecated.
1371
1372        """
1373        if sparse is None:
1374            if self.isket:
1375                sparse = False
1376            elif (self.data.nnz / (self.shape[0] * self.shape[1])) >= 0.1:
1377                sparse = False
1378        if sparse:
1379            q = Qobj()
1380            q.data, q.dims, _ = _ptrace(self, sel)
1381            return q.tidyup() if settings.auto_tidyup else q
1382        else:
1383            return _ptrace_dense(self, sel)
1384
1385    def permute(self, order):
1386        """Permutes a composite quantum object.
1387
1388        Parameters
1389        ----------
1390        order : list/array
1391            List specifying new tensor order.
1392
1393        Returns
1394        -------
1395        P : :class:`qutip.Qobj`
1396            Permuted quantum object.
1397
1398        """
1399        q = Qobj()
1400        q.data, q.dims = _permute(self, order)
1401        q.data.sort_indices()
1402        return q
1403
1404    def tidyup(self, atol=settings.auto_tidyup_atol):
1405        """Removes small elements from the quantum object.
1406
1407        Parameters
1408        ----------
1409        atol : float
1410            Absolute tolerance used by tidyup. Default is set
1411            via qutip global settings parameters.
1412
1413        Returns
1414        -------
1415        oper : :class:`qutip.Qobj`
1416            Quantum object with small elements removed.
1417
1418        """
1419        if self.data.nnz:
1420            # This does the tidyup and returns True if
1421            # The sparse data needs to be shortened
1422            if use_openmp() and self.data.nnz > 500:
1423                if omp_tidyup(self.data.data, atol, self.data.nnz,
1424                              settings.num_cpus):
1425                    self.data.eliminate_zeros()
1426            else:
1427                if cy_tidyup(self.data.data, atol, self.data.nnz):
1428                    self.data.eliminate_zeros()
1429            return self
1430        else:
1431            return self
1432
1433    def transform(self, inpt, inverse=False, sparse=True):
1434        """Basis transform defined by input array.
1435
1436        Input array can be a ``matrix`` defining the transformation,
1437        or a ``list`` of kets that defines the new basis.
1438
1439
1440        Parameters
1441        ----------
1442        inpt : array_like
1443            A ``matrix`` or ``list`` of kets defining the transformation.
1444        inverse : bool
1445            Whether to return inverse transformation.
1446        sparse : bool
1447            Use sparse matrices when possible. Can be slower.
1448
1449        Returns
1450        -------
1451        oper : :class:`qutip.Qobj`
1452            Operator in new basis.
1453
1454        Notes
1455        -----
1456        This function is still in development.
1457
1458
1459        """
1460        if isinstance(inpt, list) or (isinstance(inpt, np.ndarray) and
1461                                      len(inpt.shape) == 1):
1462            if len(inpt) != max(self.shape):
1463                raise TypeError(
1464                    'Invalid size of ket list for basis transformation')
1465            if sparse:
1466                S = sp.hstack([psi.data for psi in inpt],
1467                              format='csr', dtype=complex).conj().T
1468            else:
1469                S = np.hstack([psi.full() for psi in inpt],
1470                              dtype=complex).conj().T
1471        elif isinstance(inpt, Qobj) and inpt.isoper:
1472            S = inpt.data
1473        elif isinstance(inpt, np.ndarray):
1474            S = inpt.conj()
1475            sparse = False
1476        else:
1477            raise TypeError('Invalid operand for basis transformation')
1478
1479        # transform data
1480        if inverse:
1481            if self.isket:
1482                data = (S.conj().T) * self.data
1483            elif self.isbra:
1484                data = self.data.dot(S)
1485            else:
1486                if sparse:
1487                    data = (S.conj().T) * self.data * S
1488                else:
1489                    data = (S.conj().T).dot(self.data.dot(S))
1490        else:
1491            if self.isket:
1492                data = S * self.data
1493            elif self.isbra:
1494                data = self.data.dot(S.conj().T)
1495            else:
1496                if sparse:
1497                    data = S * self.data * (S.conj().T)
1498                else:
1499                    data = S.dot(self.data.dot(S.conj().T))
1500
1501        out = Qobj(data, dims=self.dims)
1502        out._isherm = self._isherm
1503        out.superrep = self.superrep
1504
1505        if settings.auto_tidyup:
1506            return out.tidyup()
1507        else:
1508            return out
1509
1510    def trunc_neg(self, method="clip"):
1511        """Truncates negative eigenvalues and renormalizes.
1512
1513        Returns a new Qobj by removing the negative eigenvalues
1514        of this instance, then renormalizing to obtain a valid density
1515        operator.
1516
1517
1518        Parameters
1519        ----------
1520        method : str
1521            Algorithm to use to remove negative eigenvalues. "clip"
1522            simply discards negative eigenvalues, then renormalizes.
1523            "sgs" uses the SGS algorithm (doi:10/bb76) to find the
1524            positive operator that is nearest in the Shatten 2-norm.
1525
1526        Returns
1527        -------
1528        oper : :class:`qutip.Qobj`
1529            A valid density operator.
1530
1531        """
1532        if not self.isherm:
1533            raise ValueError("Must be a Hermitian operator to remove negative "
1534                             "eigenvalues.")
1535
1536        if method not in ('clip', 'sgs'):
1537            raise ValueError("Method {} not recognized.".format(method))
1538
1539        eigvals, eigstates = self.eigenstates()
1540        if all([eigval >= 0 for eigval in eigvals]):
1541            # All positive, so just renormalize.
1542            return self.unit()
1543        idx_nonzero = eigvals != 0
1544        eigvals = eigvals[idx_nonzero]
1545        eigstates = eigstates[idx_nonzero]
1546
1547        if method == 'clip':
1548            eigvals[eigvals < 0] = 0
1549        elif method == 'sgs':
1550            eigvals = eigvals[::-1]
1551            eigstates = eigstates[::-1]
1552
1553            acc = 0.0
1554            dim = self.shape[0]
1555            n_eigs = len(eigvals)
1556
1557            for idx in reversed(range(n_eigs)):
1558                if eigvals[idx] + acc / (idx + 1) >= 0:
1559                    break
1560                else:
1561                    acc += eigvals[idx]
1562                    eigvals[idx] = 0.0
1563
1564            eigvals[:idx+1] += acc / (idx + 1)
1565
1566        return sum([
1567                val * qutip.states.ket2dm(state)
1568                for val, state in zip(eigvals, eigstates)
1569            ], Qobj(np.zeros(self.shape), dims=self.dims)
1570        ).unit()
1571
1572    def matrix_element(self, bra, ket):
1573        """Calculates a matrix element.
1574
1575        Gives the matrix element for the quantum object sandwiched between a
1576        `bra` and `ket` vector.
1577
1578        Parameters
1579        -----------
1580        bra : :class:`qutip.Qobj`
1581            Quantum object of type 'bra' or 'ket'
1582
1583        ket : :class:`qutip.Qobj`
1584            Quantum object of type 'ket'.
1585
1586        Returns
1587        -------
1588        elem : complex
1589            Complex valued matrix element.
1590
1591        Notes
1592        -----
1593        It is slightly more computationally efficient to use a ket
1594        vector for the 'bra' input.
1595        """
1596        if not self.isoper:
1597            raise TypeError("Can only get matrix elements for an operator.")
1598
1599        else:
1600            if bra.isbra and ket.isket:
1601                return zcsr_mat_elem(self.data, bra.data, ket.data, 1)
1602
1603            elif bra.isket and ket.isket:
1604                return zcsr_mat_elem(self.data, bra.data, ket.data, 0)
1605            else:
1606                err = "Can only calculate matrix elements for bra"
1607                err += " and ket vectors."
1608                raise TypeError(err)
1609
1610    def overlap(self, other):
1611        """Overlap between two state vectors or two operators.
1612
1613        Gives the overlap (inner product) between the current bra or ket Qobj
1614        and and another bra or ket Qobj. It gives the Hilbert-Schmidt overlap
1615        when one of the Qobj is an operator/density matrix.
1616
1617        Parameters
1618        -----------
1619        other : :class:`qutip.Qobj`
1620            Quantum object for a state vector of type 'ket', 'bra' or density
1621            matrix.
1622
1623        Returns
1624        -------
1625        overlap : complex
1626            Complex valued overlap.
1627
1628        Raises
1629        ------
1630        TypeError
1631            Can only calculate overlap between a bra, ket and density matrix
1632            quantum objects.
1633
1634        Notes
1635        -----
1636        Since QuTiP mainly deals with ket vectors, the most efficient inner
1637        product call is the ket-ket version that computes the product
1638        <self|other> with both vectors expressed as kets.
1639        """
1640
1641        if isinstance(other, Qobj):
1642
1643            if self.isbra:
1644                if other.isket:
1645                    return zcsr_inner(self.data, other.data, 1)
1646                elif other.isbra:
1647                    # Since we deal mainly with ket vectors, the bra-bra combo
1648                    # is not common, and not optimized.
1649                    return zcsr_inner(self.data, other.dag().data, 1)
1650                elif other.isoper:
1651                    return (qutip.states.ket2dm(self).dag() * other).tr()
1652                else:
1653                    err = "Can only calculate overlap for state vector Qobjs"
1654                    raise TypeError(err)
1655
1656            elif self.isket:
1657                if other.isbra:
1658                    return zcsr_inner(other.data, self.data, 1)
1659                elif other.isket:
1660                    return zcsr_inner(self.data, other.data, 0)
1661                elif other.isoper:
1662                    return (qutip.states.ket2dm(self).dag() * other).tr()
1663                else:
1664                    err = "Can only calculate overlap for state vector Qobjs"
1665                    raise TypeError(err)
1666
1667            elif self.isoper:
1668                if other.isket or other.isbra:
1669                    return (self.dag() * qutip.states.ket2dm(other)).tr()
1670                elif other.isoper:
1671                    return (self.dag() * other).tr()
1672                else:
1673                    err = "Can only calculate overlap for state vector Qobjs"
1674                    raise TypeError(err)
1675        raise TypeError("Can only calculate overlap for state vector Qobjs")
1676
1677    def eigenstates(self, sparse=False, sort='low', eigvals=0,
1678                    tol=0, maxiter=100000, phase_fix=None):
1679        """Eigenstates and eigenenergies.
1680
1681        Eigenstates and eigenenergies are defined for operators and
1682        superoperators only.
1683
1684        Parameters
1685        ----------
1686        sparse : bool
1687            Use sparse Eigensolver
1688
1689        sort : str
1690            Sort eigenvalues (and vectors) 'low' to high, or 'high' to low.
1691
1692        eigvals : int
1693            Number of requested eigenvalues. Default is all eigenvalues.
1694
1695        tol : float
1696            Tolerance used by sparse Eigensolver (0 = machine precision).
1697            The sparse solver may not converge if the tolerance is set too low.
1698
1699        maxiter : int
1700            Maximum number of iterations performed by sparse solver (if used).
1701
1702        phase_fix : int, None
1703            If not None, set the phase of each kets so that ket[phase_fix,0]
1704            is real positive.
1705
1706        Returns
1707        -------
1708        eigvals : array
1709            Array of eigenvalues for operator.
1710
1711        eigvecs : array
1712            Array of quantum operators representing the oprator eigenkets.
1713            Order of eigenkets is determined by order of eigenvalues.
1714
1715        Notes
1716        -----
1717        The sparse eigensolver is much slower than the dense version.
1718        Use sparse only if memory requirements demand it.
1719
1720        """
1721        evals, evecs = sp_eigs(self.data, self.isherm, sparse=sparse,
1722                               sort=sort, eigvals=eigvals, tol=tol,
1723                               maxiter=maxiter)
1724        if self.type == 'super':
1725            new_dims = [self.dims[0], [1]]
1726            new_type = 'operator-ket'
1727        else:
1728            new_dims = [self.dims[0], [1] * len(self.dims[0])]
1729            new_type = 'ket'
1730        ekets = np.empty((len(evecs),), dtype=object)
1731        ekets[:] = [Qobj(vec, dims=new_dims, type=new_type) for vec in evecs]
1732        norms = np.array([ket.norm() for ket in ekets])
1733        if phase_fix is None:
1734            phase = np.array([1] * len(ekets))
1735        else:
1736            phase = np.array([np.abs(ket[phase_fix, 0]) / ket[phase_fix, 0]
1737                              if ket[phase_fix, 0] else 1
1738                              for ket in ekets])
1739        return evals, ekets / norms * phase
1740
1741    def eigenenergies(self, sparse=False, sort='low',
1742                      eigvals=0, tol=0, maxiter=100000):
1743        """Eigenenergies of a quantum object.
1744
1745        Eigenenergies (eigenvalues) are defined for operators or superoperators
1746        only.
1747
1748        Parameters
1749        ----------
1750        sparse : bool
1751            Use sparse Eigensolver
1752        sort : str
1753            Sort eigenvalues 'low' to high, or 'high' to low.
1754        eigvals : int
1755            Number of requested eigenvalues. Default is all eigenvalues.
1756        tol : float
1757            Tolerance used by sparse Eigensolver (0=machine precision).
1758            The sparse solver may not converge if the tolerance is set too low.
1759        maxiter : int
1760            Maximum number of iterations performed by sparse solver (if used).
1761
1762        Returns
1763        -------
1764        eigvals : array
1765            Array of eigenvalues for operator.
1766
1767        Notes
1768        -----
1769        The sparse eigensolver is much slower than the dense version.
1770        Use sparse only if memory requirements demand it.
1771
1772        """
1773        return sp_eigs(self.data, self.isherm, vecs=False, sparse=sparse,
1774                       sort=sort, eigvals=eigvals, tol=tol, maxiter=maxiter)
1775
1776    def groundstate(self, sparse=False, tol=0, maxiter=100000, safe=True):
1777        """Ground state Eigenvalue and Eigenvector.
1778
1779        Defined for quantum operators or superoperators only.
1780
1781        Parameters
1782        ----------
1783        sparse : bool
1784            Use sparse Eigensolver
1785        tol : float
1786            Tolerance used by sparse Eigensolver (0 = machine precision).
1787            The sparse solver may not converge if the tolerance is set too low.
1788        maxiter : int
1789            Maximum number of iterations performed by sparse solver (if used).
1790        safe : bool (default=True)
1791            Check for degenerate ground state
1792
1793        Returns
1794        -------
1795        eigval : float
1796            Eigenvalue for the ground state of quantum operator.
1797        eigvec : :class:`qutip.Qobj`
1798            Eigenket for the ground state of quantum operator.
1799
1800        Notes
1801        -----
1802        The sparse eigensolver is much slower than the dense version.
1803        Use sparse only if memory requirements demand it.
1804
1805        """
1806        if safe:
1807            evals = 2
1808        else:
1809            evals = 1
1810        grndval, grndvec = sp_eigs(self.data, self.isherm, sparse=sparse,
1811                                   eigvals=evals, tol=tol, maxiter=maxiter)
1812        if safe:
1813            tol = 1e-15 if tol == 0 else tol
1814            if (grndval[1]-grndval[0]) <= 10*tol:
1815                print("WARNING: Ground state may be degenerate. "
1816                      "Use Q.eigenstates()")
1817        new_dims = [self.dims[0], [1] * len(self.dims[0])]
1818        grndvec = Qobj(grndvec[0], dims=new_dims)
1819        grndvec = grndvec / grndvec.norm()
1820        return grndval[0], grndvec
1821
1822    def trans(self):
1823        """Transposed operator.
1824
1825        Returns
1826        -------
1827        oper : :class:`qutip.Qobj`
1828            Transpose of input operator.
1829
1830        """
1831        out = Qobj()
1832        out.data = zcsr_transpose(self.data)
1833        out.dims = [self.dims[1], self.dims[0]]
1834        return out
1835
1836    def extract_states(self, states_inds, normalize=False):
1837        """Qobj with states in state_inds only.
1838
1839        Parameters
1840        ----------
1841        states_inds : list of integer
1842            The states that should be kept.
1843
1844        normalize : True / False
1845            Weather or not the new Qobj instance should be normalized (default
1846            is False). For Qobjs that represents density matrices or state
1847            vectors normalized should probably be set to True, but for Qobjs
1848            that represents operators in for example an Hamiltonian, normalize
1849            should be False.
1850
1851        Returns
1852        -------
1853        q : :class:`qutip.Qobj`
1854            A new instance of :class:`qutip.Qobj` that contains only the states
1855            corresponding to the indices in `state_inds`.
1856
1857        Notes
1858        -----
1859        Experimental.
1860
1861        """
1862        if self.isoper:
1863            q = Qobj(self.data[states_inds, :][:, states_inds])
1864        elif self.isket:
1865            q = Qobj(self.data[states_inds, :])
1866        elif self.isbra:
1867            q = Qobj(self.data[:, states_inds])
1868        else:
1869            raise TypeError("Can only eliminate states from operators or " +
1870                            "state vectors")
1871
1872        return q.unit() if normalize else q
1873
1874    def eliminate_states(self, states_inds, normalize=False):
1875        """Creates a new quantum object with states in state_inds eliminated.
1876
1877        Parameters
1878        ----------
1879        states_inds : list of integer
1880            The states that should be removed.
1881
1882        normalize : True / False
1883            Weather or not the new Qobj instance should be normalized (default
1884            is False). For Qobjs that represents density matrices or state
1885            vectors normalized should probably be set to True, but for Qobjs
1886            that represents operators in for example an Hamiltonian, normalize
1887            should be False.
1888
1889        Returns
1890        -------
1891        q : :class:`qutip.Qobj`
1892            A new instance of :class:`qutip.Qobj` that contains only the states
1893            corresponding to indices that are **not** in `state_inds`.
1894
1895        Notes
1896        -----
1897        Experimental.
1898
1899        """
1900        keep_indices = np.array([s not in states_inds
1901                                 for s in range(self.shape[0])]).nonzero()[0]
1902
1903        return self.extract_states(keep_indices, normalize=normalize)
1904
1905    def dnorm(self, B=None):
1906        """Calculates the diamond norm, or the diamond distance to another
1907        operator.
1908
1909        Parameters
1910        ----------
1911        B : :class:`qutip.Qobj` or None
1912            If B is not None, the diamond distance d(A, B) = dnorm(A - B)
1913            between this operator and B is returned instead of the diamond
1914            norm.
1915
1916        Returns
1917        -------
1918        d : float
1919            Either the diamond norm of this operator, or the diamond distance
1920            from this operator to B.
1921
1922        """
1923        return mts.dnorm(self, B)
1924
1925    @property
1926    def ishp(self):
1927        # FIXME: this needs to be cached in the same ways as isherm.
1928        if self.type in ["super", "oper"]:
1929            try:
1930                J = sr.to_choi(self)
1931                return J.isherm
1932            except TypeError:
1933                return False
1934        else:
1935            return False
1936
1937    @property
1938    def iscp(self):
1939        # FIXME: this needs to be cached in the same ways as isherm.
1940        if self.type in ["super", "oper"]:
1941            try:
1942                J = (
1943                    self
1944                    # We can test with either Choi or chi, since the basis
1945                    # transformation between them is unitary and hence
1946                    # preserves the CP and TP conditions.
1947                    if self.superrep in ('choi', 'chi')
1948                    else sr.to_choi(self)
1949                )
1950                # If J isn't hermitian, then that could indicate either
1951                # that J is not normal, or is normal,
1952                # but has complex eigenvalues.
1953                # In either case, it makes no sense to then demand that the
1954                # eigenvalues be non-negative.
1955                if not J.isherm:
1956                    return False
1957                eigs = J.eigenenergies()
1958                return all(eigs >= -settings.atol)
1959            except TypeError:
1960                return False
1961        else:
1962            return False
1963
1964    @property
1965    def istp(self):
1966        import qutip.superop_reps as sr
1967        if self.type in ["super", "oper"]:
1968            try:
1969                # Normalize to a super of type choi or chi.
1970                # We can test with either Choi or chi, since the basis
1971                # transformation between them is unitary and hence
1972                # preserves the CP and TP conditions.
1973                if self.type == "super" and self.superrep in ('choi', 'chi'):
1974                    qobj = self
1975                else:
1976                    qobj = sr.to_choi(self)
1977
1978                # Possibly collapse dims.
1979                if any(
1980                    len(index) > 1
1981                    for super_index in qobj.dims for index in super_index
1982                ):
1983                    qobj = Qobj(qobj, dims=collapse_dims_super(qobj.dims))
1984                else:
1985                    qobj = qobj
1986
1987                # We use the condition from John Watrous' lecture notes,
1988                # Tr_1(J(Phi)) = identity_2.
1989                tr_oper = qobj.ptrace([0])
1990                ident = ops.identity(tr_oper.shape[0])
1991                return isequal(tr_oper, ident)
1992            except TypeError:
1993                return False
1994        else:
1995            return False
1996
1997    @property
1998    def iscptp(self):
1999        from qutip.superop_reps import to_choi
2000        if self.type == "super" or self.type == "oper":
2001            reps = ('choi', 'chi')
2002            q_oper = to_choi(self) if self.superrep not in reps else self
2003            return q_oper.iscp and q_oper.istp
2004        else:
2005            return False
2006
2007    @property
2008    def isherm(self):
2009        if self._isherm is not None:
2010            # used previously computed value
2011            return self._isherm
2012
2013        self._isherm = bool(zcsr_isherm(self.data))
2014
2015        return self._isherm
2016
2017    @isherm.setter
2018    def isherm(self, isherm):
2019        self._isherm = isherm
2020
2021    def check_isunitary(self):
2022        """
2023        Checks whether qobj is a unitary matrix
2024        """
2025        if self.isoper:
2026            eye_data = fast_identity(self.shape[0])
2027            return not (
2028                np.any(
2029                    np.abs((self.data*self.dag().data - eye_data).data)
2030                    > settings.atol
2031                )
2032                or
2033                np.any(
2034                    np.abs((self.dag().data*self.data - eye_data).data)
2035                    > settings.atol
2036                )
2037            )
2038        else:
2039            return False
2040
2041    @property
2042    def isunitary(self):
2043        if self._isunitary is not None:
2044            # used previously computed value
2045            return self._isunitary
2046
2047        self._isunitary = self.check_isunitary()
2048
2049        return self._isunitary
2050
2051    @isunitary.setter
2052    def isunitary(self, isunitary):
2053        self._isunitary = isunitary
2054
2055    @property
2056    def type(self):
2057        if not self._type:
2058            self._type = type_from_dims(self.dims)
2059
2060        return self._type
2061
2062    @property
2063    def shape(self):
2064        if self.data.shape == (1, 1):
2065            return tuple([np.prod(self.dims[0]), np.prod(self.dims[1])])
2066        else:
2067            return tuple(self.data.shape)
2068
2069    @property
2070    def isbra(self):
2071        return self.type == 'bra'
2072
2073    @property
2074    def isket(self):
2075        return self.type == 'ket'
2076
2077    @property
2078    def isoperbra(self):
2079        return self.type == 'operator-bra'
2080
2081    @property
2082    def isoperket(self):
2083        return self.type == 'operator-ket'
2084
2085    @property
2086    def isoper(self):
2087        return self.type == 'oper'
2088
2089    @property
2090    def issuper(self):
2091        return self.type == 'super'
2092
2093    @staticmethod
2094    def evaluate(qobj_list, t, args):
2095        """Evaluate a time-dependent quantum object in list format. For
2096        example,
2097
2098            qobj_list = [H0, [H1, func_t]]
2099
2100        is evaluated to
2101
2102            Qobj(t) = H0 + H1 * func_t(t, args)
2103
2104        and
2105
2106            qobj_list = [H0, [H1, 'sin(w * t)']]
2107
2108        is evaluated to
2109
2110            Qobj(t) = H0 + H1 * sin(args['w'] * t)
2111
2112        Parameters
2113        ----------
2114        qobj_list : list
2115            A nested list of Qobj instances and corresponding time-dependent
2116            coefficients.
2117        t : float
2118            The time for which to evaluate the time-dependent Qobj instance.
2119        args : dictionary
2120            A dictionary with parameter values required to evaluate the
2121            time-dependent Qobj intance.
2122
2123        Returns
2124        -------
2125        output : :class:`qutip.Qobj`
2126            A Qobj instance that represents the value of qobj_list at time t.
2127
2128        """
2129
2130        q_sum = 0
2131        if isinstance(qobj_list, Qobj):
2132            q_sum = qobj_list
2133        elif isinstance(qobj_list, list):
2134            for q in qobj_list:
2135                if isinstance(q, Qobj):
2136                    q_sum += q
2137                elif (isinstance(q, list) and len(q) == 2 and
2138                      isinstance(q[0], Qobj)):
2139                    if isinstance(q[1], types.FunctionType):
2140                        q_sum += q[0] * q[1](t, args)
2141                    elif isinstance(q[1], str):
2142                        args['t'] = t
2143                        q_sum += q[0] * float(eval(q[1], globals(), args))
2144                    else:
2145                        raise TypeError('Unrecognized format for ' +
2146                                        'specification of time-dependent Qobj')
2147                else:
2148                    raise TypeError('Unrecognized format for specification ' +
2149                                    'of time-dependent Qobj')
2150        else:
2151            raise TypeError(
2152                'Unrecongized format for specification of time-dependent Qobj')
2153
2154        return q_sum
2155
2156
2157# -----------------------------------------------------------------------------
2158# This functions evaluates a time-dependent quantum object on the list-string
2159# and list-function formats that are used by the time-dependent solvers.
2160# Although not used directly in by those solvers, it can for test purposes be
2161# conventient to be able to evaluate the expressions passed to the solver for
2162# arbitrary value of time. This function provides this functionality.
2163#
2164def qobj_list_evaluate(qobj_list, t, args):
2165    """
2166    Deprecated: See Qobj.evaluate
2167    """
2168    warnings.warn("Deprecated: Use Qobj.evaluate", DeprecationWarning)
2169    return Qobj.evaluate(qobj_list, t, args)
2170
2171
2172# -----------------------------------------------------------------------------
2173#
2174# A collection of tests used to determine the type of quantum objects, and some
2175# functions for increased compatibility with quantum optics toolbox.
2176#
2177def dag(A):
2178    """Adjont operator (dagger) of a quantum object.
2179
2180    Parameters
2181    ----------
2182    A : :class:`qutip.Qobj`
2183        Input quantum object.
2184
2185    Returns
2186    -------
2187    oper : :class:`qutip.Qobj`
2188        Adjoint of input operator
2189
2190    Notes
2191    -----
2192    This function is for legacy compatibility only. It is recommended to use
2193    the ``dag()`` Qobj method.
2194
2195    """
2196    if not isinstance(A, Qobj):
2197        raise TypeError("Input is not a quantum object")
2198
2199    return A.dag()
2200
2201
2202def ptrace(Q, sel):
2203    """Partial trace of the Qobj with selected components remaining.
2204
2205    Parameters
2206    ----------
2207    Q : :class:`qutip.Qobj`
2208        Composite quantum object.
2209    sel : int/list
2210        An ``int`` or ``list`` of components to keep after partial trace.
2211
2212    Returns
2213    -------
2214    oper : :class:`qutip.Qobj`
2215        Quantum object representing partial trace with selected components
2216        remaining.
2217
2218    Notes
2219    -----
2220    This function is for legacy compatibility only. It is recommended to use
2221    the ``ptrace()`` Qobj method.
2222
2223    """
2224    if not isinstance(Q, Qobj):
2225        raise TypeError("Input is not a quantum object")
2226
2227    return Q.ptrace(sel)
2228
2229
2230def _ptrace_dense(Q, sel):
2231    rd = np.asarray(Q.dims[0], dtype=np.int32).ravel()
2232    nd = len(rd)
2233    if isinstance(sel, int):
2234        sel = np.array([sel])
2235    else:
2236        sel = np.asarray(sel)
2237    sel = list(np.sort(sel))
2238    for x in sel:
2239        if not 0 <= x < len(rd):
2240            raise IndexError("Invalid selection index in ptrace.")
2241    dkeep = (rd[sel]).tolist()
2242    qtrace = list(set(np.arange(nd)) - set(sel))
2243    dtrace = (rd[qtrace]).tolist()
2244    if len(dkeep) + len(dtrace) != len(rd):
2245        raise ValueError("Duplicate selection index in ptrace.")
2246    if not dtrace:
2247        # If we are keeping all dimensions, no need to construct an ndarray.
2248        return Q.copy()
2249    rd = list(rd)
2250    if isket(Q):
2251        vmat = (Q.full()
2252                .reshape(rd)
2253                .transpose(sel + qtrace)
2254                .reshape([np.prod(dkeep, dtype=np.int32),
2255                          np.prod(dtrace, dtype=np.int32)]))
2256        rhomat = vmat.dot(vmat.conj().T)
2257    else:
2258        rhomat = np.trace(Q.full()
2259                          .reshape(rd + rd)
2260                          .transpose(qtrace + [nd + q for q in qtrace] +
2261                                     sel + [nd + q for q in sel])
2262                          .reshape([np.prod(dtrace, dtype=np.int32),
2263                                    np.prod(dtrace, dtype=np.int32),
2264                                    np.prod(dkeep, dtype=np.int32),
2265                                    np.prod(dkeep, dtype=np.int32)]))
2266    return Qobj(rhomat, dims=[dkeep, dkeep])
2267
2268
2269def dims(inpt):
2270    """Returns the dims attribute of a quantum object.
2271
2272    Parameters
2273    ----------
2274    inpt : :class:`qutip.Qobj`
2275        Input quantum object.
2276
2277    Returns
2278    -------
2279    dims : list
2280        A ``list`` of the quantum objects dimensions.
2281
2282    Notes
2283    -----
2284    This function is for legacy compatibility only. Using the `Qobj.dims`
2285    attribute is recommended.
2286
2287    """
2288    if isinstance(inpt, Qobj):
2289        return inpt.dims
2290    else:
2291        raise TypeError("Input is not a quantum object")
2292
2293
2294def shape(inpt):
2295    """Returns the shape attribute of a quantum object.
2296
2297    Parameters
2298    ----------
2299    inpt : :class:`qutip.Qobj`
2300        Input quantum object.
2301
2302    Returns
2303    -------
2304    shape : list
2305        A ``list`` of the quantum objects shape.
2306
2307    Notes
2308    -----
2309    This function is for legacy compatibility only. Using the `Qobj.shape`
2310    attribute is recommended.
2311
2312    """
2313    if isinstance(inpt, Qobj):
2314        return Qobj.shape
2315    else:
2316        return np.shape(inpt)
2317
2318
2319def isket(Q):
2320    """
2321    Determines if given quantum object is a ket-vector.
2322
2323    Parameters
2324    ----------
2325    Q : :class:`qutip.Qobj`
2326        Quantum object
2327
2328    Returns
2329    -------
2330    isket : bool
2331        True if qobj is ket-vector, False otherwise.
2332
2333    Examples
2334    --------
2335    >>> psi = basis(5,2)
2336    >>> isket(psi)
2337    True
2338
2339    Notes
2340    -----
2341    This function is for legacy compatibility only. Using the `Qobj.isket`
2342    attribute is recommended.
2343
2344    """
2345    return True if isinstance(Q, Qobj) and Q.isket else False
2346
2347
2348def isbra(Q):
2349    """Determines if given quantum object is a bra-vector.
2350
2351    Parameters
2352    ----------
2353    Q : :class:`qutip.Qobj`
2354        Quantum object
2355
2356    Returns
2357    -------
2358    isbra : bool
2359        True if Qobj is bra-vector, False otherwise.
2360
2361    Examples
2362    --------
2363    >>> psi = basis(5,2)
2364    >>> isket(psi)
2365    False
2366
2367    Notes
2368    -----
2369    This function is for legacy compatibility only. Using the `Qobj.isbra`
2370    attribute is recommended.
2371
2372    """
2373    return True if isinstance(Q, Qobj) and Q.isbra else False
2374
2375
2376def isoperket(Q):
2377    """Determines if given quantum object is an operator in column vector form
2378    (operator-ket).
2379
2380    Parameters
2381    ----------
2382    Q : :class:`qutip.Qobj`
2383        Quantum object
2384
2385    Returns
2386    -------
2387    isoperket : bool
2388        True if Qobj is operator-ket, False otherwise.
2389
2390    Notes
2391    -----
2392    This function is for legacy compatibility only. Using the `Qobj.isoperket`
2393    attribute is recommended.
2394
2395    """
2396    return True if isinstance(Q, Qobj) and Q.isoperket else False
2397
2398
2399def isoperbra(Q):
2400    """Determines if given quantum object is an operator in row vector form
2401    (operator-bra).
2402
2403    Parameters
2404    ----------
2405    Q : :class:`qutip.Qobj`
2406        Quantum object
2407
2408    Returns
2409    -------
2410    isoperbra : bool
2411        True if Qobj is operator-bra, False otherwise.
2412
2413    Notes
2414    -----
2415    This function is for legacy compatibility only. Using the `Qobj.isoperbra`
2416    attribute is recommended.
2417
2418    """
2419    return True if isinstance(Q, Qobj) and Q.isoperbra else False
2420
2421
2422def isoper(Q):
2423    """Determines if given quantum object is a operator.
2424
2425    Parameters
2426    ----------
2427    Q : :class:`qutip.Qobj`
2428        Quantum object
2429
2430    Returns
2431    -------
2432    isoper : bool
2433        True if Qobj is operator, False otherwise.
2434
2435    Examples
2436    --------
2437    >>> a = destroy(5)
2438    >>> isoper(a)
2439    True
2440
2441    Notes
2442    -----
2443    This function is for legacy compatibility only. Using the `Qobj.isoper`
2444    attribute is recommended.
2445
2446    """
2447    return True if isinstance(Q, Qobj) and Q.isoper else False
2448
2449
2450def issuper(Q):
2451    """Determines if given quantum object is a super-operator.
2452
2453    Parameters
2454    ----------
2455    Q : :class:`qutip.Qobj`
2456        Quantum object
2457
2458    Returns
2459    -------
2460    issuper  : bool
2461        True if Qobj is superoperator, False otherwise.
2462
2463    Notes
2464    -----
2465    This function is for legacy compatibility only. Using the `Qobj.issuper`
2466    attribute is recommended.
2467
2468    """
2469    return True if isinstance(Q, Qobj) and Q.issuper else False
2470
2471
2472def isequal(A, B, tol=None):
2473    """Determines if two qobj objects are equal to within given tolerance.
2474
2475    Parameters
2476    ----------
2477    A : :class:`qutip.Qobj`
2478        Qobj one
2479    B : :class:`qutip.Qobj`
2480        Qobj two
2481    tol : float
2482        Tolerence for equality to be valid
2483
2484    Returns
2485    -------
2486    isequal : bool
2487        True if qobjs are equal, False otherwise.
2488
2489    Notes
2490    -----
2491    This function is for legacy compatibility only. Instead, it is recommended
2492    to use the equality operator of Qobj instances instead: A == B.
2493
2494    """
2495    if tol is None:
2496        tol = settings.atol
2497
2498    if not isinstance(A, Qobj) or not isinstance(B, Qobj):
2499        return False
2500
2501    if A.dims != B.dims:
2502        return False
2503
2504    Adat = A.data
2505    Bdat = B.data
2506    elems = (Adat - Bdat).data
2507    if np.any(np.abs(elems) > tol):
2508        return False
2509
2510    return True
2511
2512
2513def isherm(Q):
2514    """Determines if given operator is Hermitian.
2515
2516    Parameters
2517    ----------
2518    Q : :class:`qutip.Qobj`
2519        Quantum object
2520
2521    Returns
2522    -------
2523    isherm : bool
2524        True if operator is Hermitian, False otherwise.
2525
2526    Examples
2527    --------
2528    >>> a = destroy(4)
2529    >>> isherm(a)
2530    False
2531
2532    Notes
2533    -----
2534    This function is for legacy compatibility only. Using the `Qobj.isherm`
2535    attribute is recommended.
2536
2537    """
2538    return True if isinstance(Q, Qobj) and Q.isherm else False
2539
2540
2541# TRAILING IMPORTS
2542# We do a few imports here to avoid circular dependencies.
2543from qutip.eseries import eseries
2544import qutip.superop_reps as sr
2545import qutip.tensor as tensor
2546import qutip.operators as ops
2547import qutip.metrics as mts
2548import qutip.states
2549import qutip.superoperator
2550