class ScalarValue(Scalar):
    """Wrapper around a numeric or symbolic value

    The wrapped value may be of any of the following types::

        >>> for t in ScalarValue._val_types:
        ...     print(t)
        <class 'int'>
        <class 'float'>
        <class 'complex'>
        <class 'sympy.core.basic.Basic'>
        <class 'numpy.int64'>
        <class 'numpy.complex128'>
        <class 'numpy.float64'>

    A :class:`ScalarValue` behaves exactly like its wrapped value in all
    algebraic contexts::

        >>> 5 * ScalarValue.create(2)
        10

    Any unknown attributes or methods will be forwarded to the wrapped value
    to ensure complete "duck-typing"::

        >>> alpha = ScalarValue(sympy.symbols('alpha', positive=True))
        >>> alpha.is_positive   # same as alpha.val.is_positive
        True
        >>> ScalarValue(5).is_positive
        Traceback (most recent call last):
          ...
        AttributeError: 'int' object has no attribute 'is_positive'
    """

    @classmethod
    def create(cls, val):
        """Instatiate the :class:`ScalarValue` while recognizing :class:`Zero`
        and :class:`One`.

        :class:`Scalar` instances as `val` (including
        :class:`ScalarExpression` instances) are left unchanged. This makes
        :meth:`ScalarValue.create` a safe method for converting unknown objects
        to :class:`Scalar`.
        """
        if val in cls._invalid:
            raise ValueError("Invalid value %r" % val)
        if val == 0:
            return Zero
        elif val == 1:
            return One
        elif isinstance(val, Scalar):
            return val
        else:
            # We instantiate ScalarValue directly to avoid the overhead of
            # super().create(). Thus, there is no caching for scalars (which is
            # probably a good thing)
            return cls(val)

    def __init__(self, val):
        self._val = val
        if not isinstance(val, self._val_types):
            raise TypeError(
                "val must be one of " +
                ", ".join(["%s" % t for t in self._val_types]))
        super().__init__(val)

    def __getattr__(self, name):
        return getattr(self.val, name)

    def _diff(self, sym):
        if isinstance(self.val, sympy.Basic):
            return ScalarValue.create(sympy.diff(self.val, sym))
        else:
            return Zero

    def _simplify_scalar(self, func):
        if isinstance(self.val, sympy.Basic):
            return self.__class__.create(func(self.val))
        else:
            return self

    @property
    def val(self):
        """The wrapped scalar value"""
        return self._val

    @property
    def args(self):
        """Tuple containing the wrapped scalar value as its only element"""
        return (self._val,)

    def _series_expand(self, param, about, order):
        if isinstance(self.val, sympy.Basic):
            if about != 0:
                c = self.val.subs({param: about + param})
            else:
                c = self.val
            series = sympy.series(c, x=param, x0=0, n=None)
            res = []
            next_order = 0
            for term in series:
                c, o = term.as_coeff_exponent(param)
                if o < 0 or o.is_noninteger:
                    raise ValueError(
                        "%s is singular at expansion point %s=%s."
                        % (self, param, about))
                if o > order:
                    break
                res.extend([0] * (o - next_order))
                res.append(c)
                next_order = o + 1
            res.extend([0] * (order + 1 - next_order))
            return tuple([ScalarValue.create(c) for c in res])
        else:
            return tuple([self, ] + [Zero] * order)

    @property
    def real(self):
        """Real part"""
        if hasattr(self.val, 'real'):
            return self.val.real
        else:
            # SymPy
            return self.val.as_real_imag()[0]

    @property
    def imag(self):
        """Imaginary part"""
        if hasattr(self.val, 'imag'):
            return self.val.imag
        else:
            # SymPy
            return self.val.as_real_imag()[1]

    def _adjoint(self):
        return self.__class__(self.val.conjugate())

    def __eq__(self, other):
        if isinstance(other, ScalarValue):
            return self.val == other.val
        else:
            return self.val == other

    def __lt__(self, other):
        if isinstance(other, ScalarValue):
            return self.val < other.val
        else:
            return self.val < other

    def __le__(self, other):
        if isinstance(other, ScalarValue):
            return self.val <= other.val
        else:
            return self.val <= other

    def __gt__(self, other):
        if isinstance(other, ScalarValue):
            return self.val > other.val
        else:
            return self.val > other

    def __ge__(self, other):
        if isinstance(other, ScalarValue):
            return self.val >= other.val
        else:
            return self.val >= other

    def __hash__(self):
        return hash(self.val)

    def __neg__(self):
        return self.create(-self.val)

    def __abs__(self):
        return self.create(abs(self.val))

    def __add__(self, other):
        if isinstance(other, ScalarValue):
            return self.create(self.val + other.val)
        elif isinstance(other, self._val_types):
            return self.create(self.val + other)
        elif other == 1:
            return self.create(self.val + 1)
        else:
            return super().__add__(other)  # other == 0

    def __sub__(self, other):
        if isinstance(other, ScalarValue):
            return self.create(self.val - other.val)
        elif isinstance(other, self._val_types):
            return self.create(self.val - other)
        elif other == 1:
            return self.create(self.val - 1)
        else:
            return super().__sub__(other)  # other == 0

    def __mul__(self, other):
        if isinstance(other, ScalarValue):
            return self.create(self.val * other.val)
        elif isinstance(other, self._val_types):
            return self.create(self.val * other)
        else:
            return super().__mul__(other)  # other == 0, 1

    def __floordiv__(self, other):
        if isinstance(other, ScalarValue):
            return self.create(self.val // other.val)
        elif isinstance(other, self._val_types):
            return self.create(self.val // other)
        else:
            return super().__floordiv__(other)  # other == 0, 1

    def __truediv__(self, other):
        if isinstance(other, ScalarValue):
            return self.create(self.val / other.val)
        elif isinstance(other, self._val_types):
            try:
                return self.create(self.val / other)
            except ValueError:
                # sympy may produce 'infinity', which `create` catches as a
                # ValueError
                raise ZeroDivisionError("integer division or modulo by zero")
        else:
            return super().__truediv__(other)  # other == 0, 1

    def __mod__(self, other):
        if isinstance(other, ScalarValue):
            return self.create(self.val % other.val)
        elif isinstance(other, self._val_types):
            return self.create(self.val % other)
        else:
            return super().__mod__(other)  # other == 0, 1

    def __pow__(self, other):
        if isinstance(other, ScalarValue):
            return self.create(self.val**other.val)
        elif isinstance(other, self._val_types):
            return self.create(self.val**other)
        else:
            return super().__pow__(other)  # other == 0, 1

    def __radd__(self, other):
        if other == 1:
            return self.create(1 + self.val)
        elif isinstance(other, self._val_types):
            return self.create(other + self.val)
        else:
            return super().__radd__(other)  # other == 0

    def __rsub__(self, other):
        if other == 1:
            return self.create(1 - self.val)
        elif isinstance(other, self._val_types):
            return self.create(other - self.val)
        else:
            return super().__radd__(other)  # other == 0

    def __rmul__(self, other):
        if isinstance(other, self._val_types):
            return self.create(other * self.val)
        else:
            return super().__rmul__(other)  # other == 0, 1

    def __rfloordiv__(self, other):
        if other == 1:
            return self.create(1 // self.val)
        elif isinstance(other, self._val_types):
            return self.create(other // self.val)
        else:
            return super().__rfloordiv__(other)  # other == 0

    def __rtruediv__(self, other):
        if other == 1:
            return self.create(1 / self.val)
        elif isinstance(other, self._val_types):
            return self.create(other / self.val)
        else:
            return super().__rtruediv__(other)  # other == 0, 1/x -> x^(-1)

    def __rmod__(self, other):
        if isinstance(other, self._val_types):
            return self.create(other % self.val)
        else:
            return super().__rmod__(other)  # other == 0

    def __rpow__(self, other):
        if isinstance(other, self._val_types):
            return self.create(other**self.val)
        else:
            return super().__rpow__(other)  # other == 0, 1

    def __complex__(self):
        return complex(self.val)

    def __int__(self):
        return int(self.val)

    def __float__(self):
        return float(self.val)

    def _sympy_(self):
        return sympy.sympify(self.val) def is_scalar(scalar):
    """Check if `scalar` is a :class:`Scalar` or a scalar value

    Specifically, whether `scalar` is an instance of :class:`Scalar` or an
    instance of a numeric or symbolic type that could be wrapped in
    :class:`ScalarValue`.

    For internal use only.
    """
    return isinstance(scalar, Scalar) or isinstance(scalar, Scalar._val_types) class AlgebraError(AlgebraException):
    """Base class for all algebraic errors"""
    pass class SLH(Circuit, Expression):
    """Element of the SLH algebra

    The SLH class encapsulate an open system model that is parametrized the a
    scattering matrix (S), a column vector of Lindblad operators (L), and a
    Hamiltonian (H).

    Args:
        S (Matrix): The scattering matrix (with in general Operator-valued
            elements)
        L (Matrix): The coupling vector (with in general Operator-valued
            elements)
        H (Operator): The internal Hamiltonian operator
    """

    def __init__(self, S, L, H):
        if not isinstance(S, Matrix):
            S = Matrix(S)
        if not isinstance(L, Matrix):
            L = Matrix(L)
        if S.shape[0] != L.shape[0]:
            raise ValueError('S and L misaligned: S = {!r}, L = {!r}'
                             .format(S, L))
        if L.shape[1] != 1:
            raise ValueError(("L has wrong shape %s. L must be a column vector "
                              "of operators (shape n × 1)") % str(L.shape))

        if not all(isinstance(s, Operator) for s in S.matrix.ravel()):
            S = S * IdentityOperator
        if not all(isinstance(l, Operator) for l in L.matrix.ravel()):
            L = L * IdentityOperator
        if not isinstance(H, Operator):
            H = H * IdentityOperator

        self.S = S  #: Scattering matrix
        self.L = L  #: Coupling vector
        self.H = H  #: Hamiltonian
        super().__init__(S, L, H)

    @property
    def args(self):
        return self.S, self.L, self.H

    @property
    def Ls(self):
        """Lindblad operators (entries of the L vector), as a list"""
        return list(self.L.matrix[:, 0])

    @property
    def cdim(self):
        """The circuit dimension"""
        return self.S.shape[0]

    def _creduce(self):
        return self

    @property
    def space(self):
        """Total Hilbert space"""
        args_spaces = (self.S.space, self.L.space, self.H.space)
        return ProductSpace.create(*args_spaces)

    @property
    def free_symbols(self):
        """Set of all symbols occcuring in S, L, or H"""
        return set.union(
            self.S.free_symbols, self.L.free_symbols, self.H.free_symbols)

    def series_with_slh(self, other):
        """Series product with another :class:`SLH` object

        Args:
            other (SLH): An upstream SLH circuit.

        Returns:
            SLH: The combined system.
        """
        new_S = self.S * other.S
        new_L = self.S * other.L + self.L

        def ImAdjoint(m):
            return (m.H - m) * (I / 2)

        delta = ImAdjoint(self.L.adjoint() * self.S * other.L)

        if isinstance(delta, Matrix):
            new_H = self.H + other.H + delta[0, 0]
        else:
            assert delta == 0
            new_H = self.H + other.H

        return SLH(new_S, new_L, new_H)

    def concatenate_slh(self, other):
        """Concatenation with another :class:`SLH` object"""
        selfS = self.S
        otherS = other.S
        new_S = block_matrix(
                selfS, zerosm((selfS.shape[0], otherS.shape[1]), dtype=int),
                zerosm((otherS.shape[0], selfS.shape[1]), dtype=int), otherS)
        new_L = vstackm((self.L, other.L))
        new_H = self.H + other.H

        return SLH(new_S, new_L, new_H)

    def _toSLH(self):
        return self

    def expand(self):
        """Expand out all operator expressions within S, L and H

        Return a new :class:`SLH` object with these expanded expressions.
        """
        return SLH(self.S.expand(), self.L.expand(), self.H.expand())

    def simplify_scalar(self, func=sympy.simplify):
        """Simplify all scalar expressions within S, L and H

        Return a new :class:`SLH` object with the simplified expressions.

        See also: :meth:`.QuantumExpression.simplify_scalar`
        """
        return SLH(
            self.S.simplify_scalar(func=func),
            self.L.simplify_scalar(func=func),
            self.H.simplify_scalar(func=func))

    def _series_inverse(self):
        return SLH(self.S.adjoint(), - self.S.adjoint() * self.L, -self.H)

    def _feedback(self, *, out_port, in_port):

        if not isinstance(self.S, Matrix) or not isinstance(self.L, Matrix):
            return Feedback(self, out_port=out_port, in_port=in_port)

        from sympy.core.numbers import ComplexInfinity, Infinity
        sympyOne = sympify(1)

        n = self.cdim - 1

        if out_port != n:
            return (
                map_channels({out_port: n}, self.cdim).toSLH() <<
                self
            ).feedback(in_port=in_port)
        elif in_port != n:
            return (
                self << map_channels({n: in_port}, self.cdim).toSLH()
            ).feedback()

        S, L, H = self.S, self.L, self.H

        one_minus_Snn = sympyOne - S[n, n]

        if isinstance(one_minus_Snn, Operator):
            if one_minus_Snn is IdentityOperator:
                one_minus_Snn = 1
            elif (isinstance(one_minus_Snn, ScalarTimesOperator) and
                  one_minus_Snn.term is IdentityOperator):
                one_minus_Snn = one_minus_Snn.coeff
            else:
                raise AlgebraError('Inversion not implemented for general'
                                   ' operators: {}'.format(one_minus_Snn))

        one_minus_Snn_inv = sympyOne / one_minus_Snn
        if one_minus_Snn_inv in [Infinity, ComplexInfinity]:
            raise AlgebraError(
                "Ill-posed network: singularity in feedback [%s]%d->%d"
                % (str(self), out_port, in_port))

        new_S = S[:n, :n] + S[:n, n:] * one_minus_Snn_inv * S[n:, :n]
        new_L = L[:n] + S[:n, n] * one_minus_Snn_inv * L[n]

        def ImAdjoint(m):
            return (m.H - m) * (I / 2)

        delta_H = ImAdjoint(
            (L.adjoint() * S[:, n:]) * one_minus_Snn_inv * L[n, 0])

        if isinstance(delta_H, Matrix):
            delta_H = delta_H[0, 0]
        new_H = H + delta_H

        return SLH(new_S, new_L, new_H)

    def symbolic_liouvillian(self):
        from qnet.algebra.core.super_operator_algebra import liouvillian

        return liouvillian(self.H, self.L)

    def symbolic_master_equation(self, rho=None):
        """Compute the symbolic Liouvillian acting on a state rho

        If no rho is given, an OperatorSymbol is created in its place.
        This correspnds to the RHS of the master equation
        in which an average is taken over the external noise degrees of
        freedom.

        Args:
            rho (Operator): A symbolic density matrix operator

        Returns:
            Operator: The RHS of the master equation.

        """
        L, H = self.L, self.H
        if rho is None:
            rho = OperatorSymbol('rho', hs=self.space)
        return (-I * (H * rho - rho * H) +
                sum(Lk * rho * adjoint(Lk) -
                    (adjoint(Lk) * Lk * rho + rho * adjoint(Lk) * Lk) / 2
                    for Lk in L.matrix.ravel()))

    def symbolic_heisenberg_eom(
            self, X=None, noises=None, expand_simplify=True):
        """Compute the symbolic Heisenberg equations of motion of a system
        operator X.  If no X is given, an OperatorSymbol is created in its
        place.  If no noises are given, this correspnds to the
        ensemble-averaged Heisenberg equation of motion.

        Args:
            X (Operator): A system operator
            noises (Operator): A vector of noise inputs

        Returns:
            Operator: The RHS of the Heisenberg equations of motion of X.
        """
        L, H = self.L, self.H

        if X is None:
            X = OperatorSymbol('X', hs=(L.space | H.space))

        summands = [I * (H * X - X * H), ]
        for Lk in L.matrix.ravel():
            summands.append(adjoint(Lk) * X * Lk)
            summands.append(-(adjoint(Lk) * Lk * X + X * adjoint(Lk) * Lk) / 2)

        if noises is not None:
            if not isinstance(noises, Matrix):
                noises = Matrix(noises)
            LambdaT = (noises.adjoint().transpose() * noises.transpose()).transpose()
            assert noises.shape == L.shape
            S = self.S
            summands.append((adjoint(noises) * S.adjoint() * (X * L - L * X))
                            .expand()[0, 0])
            summand = (((L.adjoint() * X - X * L.adjoint()) * S * noises)
                       .expand()[0, 0])
            summands.append(summand)
            if len(S.space & X.space):
                comm = (S.adjoint() * X * S - X)
                summands.append((comm * LambdaT).expand().trace())

        ret = OperatorPlus.create(*summands)
        if expand_simplify:
            ret = ret.expand().simplify_scalar()
        return ret

    def __iter__(self):
        return iter((self.S, self.L, self.H))

    def _coherent_input(self, *input_amps):
        return super(SLH, self)._coherent_input(*input_amps).toSLH() def move_drive_to_H(slh, which=None, expand_simplify=True):
    r'''Move coherent drives from the Lindblad operators to the Hamiltonian.

    For the given SLH model, move inhomogeneities in the Lindblad operators (resulting
    from the presence of a coherent drive, see :class:`CoherentDriveCC`) to the
    Hamiltonian.

    This exploits the invariance of the Lindblad master equation under the
    transformation  (cf. Breuer and Pettrucione, Ch 3.2.1)

    .. math::
        :nowrap:

        \begin{align}
            \Op{L}_i &\longrightarrow \Op{L}_i' = \Op{L}_i - \alpha_i  \\
            \Op{H}   &\longrightarrow
            \Op{H}' = \Op{H} + \frac{1}{2i} \sum_j
                    (\alpha_j \Op{L}_j^{\dagger} - \alpha_j^* \Op{L}_j)
        \end{align}

    In the context of SLH, this transformation is achieved by feeding `slh` into

    .. math::
        \SLH(\identity, -\mat{\alpha}, 0)

    where $\mat{\alpha}$ has the elements $\alpha_i$.

    Parameters
    ----------
    slh : SLH
        SLH model to transform. If `slh` does not contain any inhomogeneities, it is
        invariant under the transformation.

    which : sequence or None
        Sequence of circuit dimensions to apply the transform to. If None, all
        dimensions are transformed.

    expand_simplify : bool
        if True, expand and simplify the new SLH object before returning. This has no
        effect if `slh` does not contain any inhomogeneities.

    Returns
    -------
    new_slh : SLH
        Transformed SLH model.

    '''
    if which is None:
        which = []
    scalarcs = []
    for jj, L in enumerate(slh.Ls):
        if not which or jj in which:
            scalarcs.append(-get_coeffs(L.expand())[IdentityOperator])
        else:
            scalarcs.append(0)

    if np.all(np.array(scalarcs) == 0):
        return slh
    new_slh = SLH(identity_matrix(slh.cdim), scalarcs, 0) << slh
    if expand_simplify:
        return new_slh.expand().simplify_scalar()
    return new_slh class Operation(Expression, metaclass=ABCMeta):
    """Base class for "operations"

    Operations are Expressions that act algebraically on other expressions
    (their "operands").

    Operations differ from more general Expressions by the convention that the
    arguments of the Operator are exactly the operands (which must be members
    of the algebra!) Any other parameters (non-operands) that may be required
    must be given as keyword-arguments.
    """

    def __init__(self, *operands, **kwargs):
        self._operands = operands
        super().__init__(*operands, **kwargs)

    @property
    def operands(self):
        """Tuple of operands of the operation"""
        return self._operands

    @property
    def args(self):
        """Alias for operands"""
        return self._operands class Operator(QuantumExpression, metaclass=ABCMeta):
    """Base class for all quantum operators."""

    def pseudo_inverse(self):
        """Pseudo-inverse $\Op{X}^+$ of the operator $\Op{X}$

        It is defined via the relationship

        .. math::

            \Op{X} \Op{X}^+ \Op{X} =  \Op{X}  \\
            \Op{X}^+ \Op{X} \Op{X}^+ =  \Op{X}^+  \\
            (\Op{X}^+ \Op{X})^\dagger = \Op{X}^+ \Op{X}  \\
            (\Op{X} \Op{X}^+)^\dagger = \Op{X} \Op{X}^+
        """
        return self._pseudo_inverse()

    @abstractmethod
    def _pseudo_inverse(self):
        raise NotImplementedError(self.__class__.__name__)

    def expand_in_basis(self, basis_states=None, hermitian=False):
        """Write the operator as an expansion into all
        :class:`KetBras <.KetBra>`
        spanned by `basis_states`.

        Args:
            basis_states (list or None): List of basis states (:class:`.State`
                instances) into which to expand the operator. If None, use the
                operator's `space.basis_states`
            hermitian (bool): If True, assume that the operator is Hermitian
                and represent all elements in the lower triangle of the
                expansion via :class:`OperatorPlusMinusCC`. This is meant to
                enhance readability

        Raises:
            .BasisNotSetError: If `basis_states` is None and the operator's
                Hilbert space has no well-defined basis

        Example:

            >>> hs = LocalSpace(1, basis=('g', 'e'))
            >>> op = LocalSigma('g', 'e', hs=hs) + LocalSigma('e', 'g', hs=hs)
            >>> print(ascii(op, sig_as_ketbra=False))
            sigma_e,g^(1) + sigma_g,e^(1)
            >>> print(ascii(op.expand_in_basis()))
            |e><g|^(1) + |g><e|^(1)
            >>> print(ascii(op.expand_in_basis(hermitian=True)))
            |g><e|^(1) + c.c.
        """
        from qnet.algebra.core.state_algebra import KetBra
        # KetBra is imported locally to avoid circular imports
        if basis_states is None:
            basis_states = list(self.space.basis_states)
        else:
            basis_states = list(basis_states)
        diag_terms = []
        terms = []
        for i, ket_i in enumerate(basis_states):
            for j, ket_j in enumerate(basis_states):
                if i > j and hermitian:
                    continue
                op_ij = (ket_i.dag() * self * ket_j).expand()
                ketbra = KetBra(ket_i, ket_j)
                term = op_ij * ketbra
                if term is not ZeroOperator:
                    if i == j:
                        diag_terms.append(op_ij * ketbra)
                    else:
                        terms.append(op_ij * ketbra)
        if hermitian:
            res = OperatorPlus.create(*diag_terms)
            if len(terms) > 0:
                res = res + OperatorPlusMinusCC(OperatorPlus.create(*terms))
            return res
        else:
            return (
                OperatorPlus.create(*diag_terms) +
                OperatorPlus.create(*terms)) class IdentityOperator(Operator, metaclass=Singleton):
    """``IdentityOperator`` constant (singleton) object."""

    _order_index = 2

    @property
    def space(self):
        """:class:`.TrivialSpace`"""
        return TrivialSpace

    @property
    def args(self):
        return tuple()

    def _diff(self, sym):
        return ZeroOperator

    def _adjoint(self):
        return self

    def _pseudo_inverse(self):
        return self class ZeroOperator(Operator, metaclass=Singleton):
    """``ZeroOperator`` constant (singleton) object."""

    _order_index = 2

    @property
    def space(self):
        """:class:`.TrivialSpace`"""
        return TrivialSpace

    @property
    def args(self):
        return tuple()

    def _diff(self, sym):
        return self

    def _adjoint(self):
        return self

    def _pseudo_inverse(self):
        return self class LocalOperator(Operator, metaclass=ABCMeta):
    """Base class for "known" operators on a :class:`LocalSpace`

    All :class:`LocalOperator` instances have known algebraic properties and a
    fixed associated identifier (symbol) that is used when printing that
    operator. A custom identifier can be used through the associated
    :class:`.LocalSpace`'s `local_identifiers` parameter. For example::

        >>> hs1_custom = LocalSpace(1, local_identifiers={'Destroy': 'b'})
        >>> b = Destroy(hs=hs1_custom)
        >>> ascii(b)
        'b^(1)'

    Note:
        It is recommended that subclasses use the :func:`.properties_for_args`
        class decorator if they define any position arguments (via the
        :attr:`_arg_names` class attribute)
    """

    simplifications = [implied_local_space(keys=['hs', ]), ]

    _identifier = None  # must be overridden by subclasses!
    _dagger = False  # do representations include a dagger?
    _arg_names = ()  # names of args that can be passed to __init__
    _scalar_args = True  # convert args to Scalar?
    _hs_cls = LocalSpace  # allowed type of `hs`
    _rx_identifier = re.compile('^[A-Za-z][A-Za-z0-9]*(_[A-Za-z0-9().+-]+)?$')

    def __init__(self, *args, hs):
        if len(args) != len(self._arg_names):
            raise ValueError("expected %d arguments, gotten %d"
                             % (len(self._arg_names), len(args)))
        if self._scalar_args:
            args = [ScalarValue.create(arg) for arg in args]
        for i, arg_name in enumerate(self._arg_names):
            self.__dict__['_%s' % arg_name] = args[i]
        hs = ensure_local_space(hs, cls=self._hs_cls)
        self._hs = hs
        if self._identifier is None:
            raise TypeError(
                r"Can't instantiate abstract class %s with undefined "
                r"_identifier" % self.__class__.__name__)
        self._args = args
        super().__init__(*args, hs=hs)

    @property
    def space(self):
        """Hilbert space of the operator (:class:`.LocalSpace` instance)"""
        return self._hs

    @property
    def args(self):
        """The positional arguments used for instantiating the operator"""
        return tuple(self._args)

    @property
    def kwargs(self):
        """The keyword arguments used for instantiating the operator"""
        return OrderedDict([('hs', self._hs)])

    @property
    def identifier(self):
        """The identifier (symbol) that is used when printing the operator.

        A custom identifier can be used through the associated
        :class:`.LocalSpace`'s `local_identifiers` parameter. For example::

            >>> a = Destroy(hs=1)
            >>> a.identifier
            'a'
            >>> hs1_custom = LocalSpace(1, local_identifiers={'Destroy': 'b'})
            >>> b = Destroy(hs=hs1_custom)
            >>> b.identifier
            'b'
            >>> ascii(b)
            'b^(1)'
        """

        identifier = self._hs._local_identifiers.get(
            self.__class__.__name__, self._identifier)
        if not self._rx_identifier.match(identifier):
            raise ValueError(
                "identifier '%s' does not match pattern '%s'"
                % (identifier, self._rx_identifier.pattern))
        return identifier

    def _diff(self, sym):
        return OperatorDerivative(self, derivs={sym: 1})

    def _simplify_scalar(self, func):
        if self._scalar_args:
            args = [arg.simplify_scalar(func=func) for arg in self.args]
            return self.create(*args, hs=self.space)
        else:
            return super()._simplify_scalar(func=func) class LocalSigma(LocalOperator):
    r'''Level flip operator between two levels of a :class:`.LocalSpace`

    .. math::

        \Op{\sigma}_{jk}^{\rm hs} =
        \left| j\right\rangle_{\rm hs} \left \langle k \right |_{\rm hs}

    For $j=k$ this becomes a projector $\Op{P}_k$ onto the eigenstate
    $\ket{k}$; see :class:`LocalProjector`.

    Args:
        j (int or str): The label or index identifying $\ket{j}$
        k (int or str):  The label or index identifying $\ket{k}$
        hs (LocalSpace or int or str): The Hilbert space on which the
            operator acts. If an :class:`int` or a :class:`str`, an implicit
            Hilbert space will be constructed as a subclass of
            :class:`LocalSpace`, as configured by :func:`init_algebra`.

    Note:
        The parameters `j` or `k` may be an integer or a string. A string
        refers to the label of an eigenstate in the basis of `hs`, which needs
        to be set. An integer refers to the (zero-based) index of eigenstate of
        the Hilbert space. This works if `hs` has an unknown dimension.
        Assuming the Hilbert space has a defined basis, using integer or string
        labels is equivalent::

            >>> hs = LocalSpace('tls', basis=('g', 'e'))
            >>> LocalSigma(0, 1, hs=hs) == LocalSigma('g', 'e', hs=hs)
            True

    Raises:
        ValueError: If `j` or `k` are invalid value for the given `hs`

    Printers should represent this operator either in braket notation, or using
    the operator identifier

        >>> LocalSigma(0, 1, hs=0).identifier
        'sigma'

    For ``j == k``, an alternative (fixed) identifier may be used

        >>> LocalSigma(0, 0, hs=0)._identifier_projector
        'Pi'
    '''

    _identifier = "sigma"
    _identifier_projector = "Pi"
    _rx_identifier = re.compile('^[A-Za-z][A-Za-z0-9]*$')
    _arg_names = ('j', 'k')
    _scalar_args = False  # args are labels, not scalar coefficients
    _rules = OrderedDict()
    simplifications = [implied_local_space(keys=['hs', ]), match_replace]

    def __init__(self, j, k, *, hs):
        if isinstance(hs, (str, int)):
            hs = self._default_hs_cls(hs)
        hs._unpack_basis_label_or_index(j)  # for applying checks only ...
        hs._unpack_basis_label_or_index(k)  # ... (disregard returned tuple)
        if hs.has_basis:
            # normalize integer i/j to str label, if possible
            if isinstance(j, int):
                j = hs.basis_labels[j]
            if isinstance(k, int):
                k = hs.basis_labels[k]
        super().__init__(j, k, hs=hs)

    @property
    def args(self):
        """The two eigenstate labels `j` and `k` that the operator connects"""
        return self.j, self.k

    @property
    def index_j(self):
        """Index `j` or (zero-based) index of the label `j` in the basis"""
        if isinstance(self.j, (int, SymbolicLabelBase)):
            return self.j
        else:
            try:
                return self.space.basis_labels.index(self.j)
            except ValueError:
                raise ValueError(
                    "%r is not one of the basis labels %r"
                    % (self.j, self.space.basis_labels))

    @property
    def index_k(self):
        """Index `k` or (zero-based) index of the label `k` in the basis"""
        if isinstance(self.k, (int, SymbolicLabelBase)):
            return self.k
        else:
            try:
                return self.space.basis_labels.index(self.k)
            except ValueError:
                raise ValueError(
                    "%r is not one of the basis labels %r"
                    % (self.k, self.space.basis_labels))

    def raise_jk(self, j_incr=0, k_incr=0):
        r'''Return a new :class:`LocalSigma` instance with incremented `j`,
        `k`, on the same Hilbert space:

        .. math::

            \Op{\sigma}_{jk}^{\rm hs} \rightarrow \Op{\sigma}_{j'k'}^{\rm hs}

        This is the result of multiplying $\Op{\sigma}_{jk}^{\rm hs}$
        with any raising or lowering operators.

        If $j'$ or $k'$ are outside the Hilbert space ${\rm hs}$, the result is
        the :obj:`ZeroOperator` .

        Args:
            j_incr (int): The increment between labels $j$ and $j'$
            k_incr (int): The increment between labels $k$ and $k'$. Both
                increments may be negative.
        '''
        try:
            if isinstance(self.j, int):
                new_j = self.j + j_incr
            else:  # str
                new_j = self.space.next_basis_label_or_index(self.j, j_incr)
            if isinstance(self.k, int):
                new_k = self.k + k_incr
            else:  # str or SymbolicLabelBase
                new_k = self.space.next_basis_label_or_index(self.k, k_incr)
            return LocalSigma.create(new_j, new_k, hs=self.space)
        except (IndexError, ValueError):
            return ZeroOperator

    def _adjoint(self):
        return LocalSigma(j=self.k, k=self.j, hs=self.space)

    def _pseudo_inverse(self):
        return self._adjoint() class OperatorPlus(QuantumPlus, Operator):
    """Sum of Operators"""

    _neutral_element = ZeroOperator
    _binary_rules = OrderedDict()
    simplifications = [
        assoc, scalars_to_op, orderby, collect_summands, match_replace_binary]

    def _pseudo_inverse(self):
        return PseudoInverse(self) class OperatorTimes(QuantumTimes, Operator):
    """Product of operators

    This serves both as a product within a Hilbert space as well as a tensor
    product."""

    _neutral_element = IdentityOperator
    _binary_rules = OrderedDict()
    simplifications = [assoc, orderby, filter_neutral, match_replace_binary]

    def _pseudo_inverse(self):
        return self.__class__.create(
                *[o._pseudo_inverse() for o in reversed(self.operands)]) class ScalarTimesOperator(Operator, ScalarTimesQuantumExpression):
    """Product of a :class:`.Scalar` coefficient and an :class:`Operator`"""

    _rules = OrderedDict()
    simplifications = [match_replace, ]

    @staticmethod
    def has_minus_prefactor(c):
        """
        For a scalar object c, determine whether it is prepended by a "-" sign.
        """
        # TODO: check if this is necessary; if yes, move
        cs = str(c).strip()
        return cs[0] == "-"

    def _pseudo_inverse(self):
        c, t = self.operands
        return t.pseudo_inverse() / c

    def __eq__(self, other):
        # TODO: review, and add this to ScalarTimesQuantumExpression
        if self.term is IdentityOperator and is_scalar(other):
            return self.coeff == other
        return super().__eq__(other)

    def __hash__(self):
        # TODO: review, and add this to ScalarTimesQuantumExpression
        return super().__hash__()

    def _adjoint(self):
        return ScalarTimesOperator(self.coeff.conjugate(), self.term.adjoint()) class Adjoint(QuantumAdjoint, Operator):
    """Symbolic Adjoint of an operator"""

    simplifications = [
        scalars_to_op, delegate_to_method('_adjoint')]
    # The reason that Adjoint does not have have `match_replace` in
    # `simplifications`, respectively a `_rules` class attribute is that the
    # `_adjoint` property that we delegate to is mandatory. Thus, if we had
    # rules on top of that, it would create the confusing situation of the rule
    # contradicting the `_adjoint` property.

    def _pseudo_inverse(self):
        return self.operand.pseudo_inverse().adjoint() class PseudoInverse(SingleQuantumOperation, Operator):
    r"""Unevaluated pseudo-inverse $\Op{X}^+$ of an operator $\Op{X}$

    It is defined via the relationship

    .. math::

        \Op{X} \Op{X}^+ \Op{X} =  \Op{X}  \\
        \Op{X}^+ \Op{X} \Op{X}^+ =  \Op{X}^+  \\
        (\Op{X}^+ \Op{X})^\dagger = \Op{X}^+ \Op{X}  \\
        (\Op{X} \Op{X}^+)^\dagger = \Op{X} \Op{X}^+
    """
    simplifications = [
        scalars_to_op, delegate_to_method('_pseudo_inverse')]
    # `PseudoInverse` does not use rules because it delegates to
    # `_pseudo_inverse`, cf. `Adjoint`

    def _expand(self):
        return self

    def _pseudo_inverse(self):
        return self.operand

    def _adjoint(self):
        return Adjoint(self) class OperatorTrace(SingleQuantumOperation, Operator):
    r'''(Partial) trace of an operator

    Trace of an operator `op` ($\Op{O}) over the degrees
    of freedom of a Hilbert space `over_space` ($\mathcal{H}$):

    .. math::

        {\rm Tr}_{\mathcal{H}} \Op{O}

    Args:
        over_space (.HilbertSpace): The degrees of freedom to trace over
        op (Operator): The operator to take the trace of.
    '''
    _rules = OrderedDict()
    simplifications = [
        scalars_to_op, implied_local_space(keys=['over_space', ]),
        match_replace, ]

    def __init__(self, op, *, over_space):
        if isinstance(over_space, (int, str)):
            over_space = self._default_hs_cls(over_space)
        assert isinstance(over_space, HilbertSpace)
        self._over_space = over_space
        super().__init__(op, over_space=over_space)
        self._space = None

    @property
    def kwargs(self):
        return {'over_space': self._over_space}

    @property
    def operand(self):
        return self.operands[0]

    @property
    def space(self):
        if self._space is None:
            return self.operands[0].space / self._over_space
        return self._space

    def _expand(self):
        s = self._over_space
        o = self.operand
        return OperatorTrace.create(o.expand(), over_space=s)

    def _series_expand(self, param, about, order):
        ope = self.operand.series_expand(param, about, order)
        return tuple(OperatorTrace.create(opet, over_space=self._over_space)
                     for opet in ope)

    def _diff(self, sym):
        s = self._over_space
        o = self.operand
        return OperatorTrace.create(o._diff(sym), over_space=s)

    def _adjoint(self):
        # there is a rule Tr[A^\dagger] -> Tr[A]^\dagger, which we don't want
        # to counteract here with an inverse rule
        return Adjoint(self)

    def _pseudo_inverse(self):
        return PseudoInverse(self) class NullSpaceProjector(SingleQuantumOperation, Operator):
    r"""Projection operator onto the nullspace of its operand

    Returns the operator :math:`\mathcal{P}_{{\rm Ker} X}` with

    .. math::

        X \mathcal{P}_{{\rm Ker} X}
          = 0
        \Leftrightarrow
        X (1 - \mathcal{P}_{{\rm Ker} X})
          = X \\
        \mathcal{P}_{{\rm Ker} X}^\dagger
          = \mathcal{P}_{{\rm Ker} X}
          = \mathcal{P}_{{\rm Ker} X}^2
    """

    _rules = OrderedDict()
    simplifications = [scalars_to_op, match_replace, ]

    def _expand(self):
        return self

    def _adjoint(self):
        return Adjoint(self)

    def _pseudo_inverse(self):
        return PseudoInverse(self) class Jz(SpinOperator):
    """Spin (angular momentum) operator in z-direction

    $\Op{J}_z$ is the $z$ component of a general spin operator acting
    on a particular :class:`SpinSpace` `hs` of freedom with well defined spin
    quantum number $s$.  It is Hermitian::

        >>> hs = SpinSpace(1, spin=(1, 2))
        >>> print(ascii(Jz(hs=hs).adjoint()))
        J_z^(1)

    :class:`Jz`, :class:`Jplus` and :class:`Jminus` satisfy the angular
    momentum commutator algebra::

        >>> print(ascii((Jz(hs=hs) * Jplus(hs=hs) -
        ...              Jplus(hs=hs)*Jz(hs=hs)).expand()))
        J_+^(1)

        >>> print(ascii((Jz(hs=hs) * Jminus(hs=hs) -
        ...              Jminus(hs=hs)*Jz(hs=hs)).expand()))
        -J_-^(1)

        >>> print(ascii((Jplus(hs=hs) * Jminus(hs=hs)
        ...              - Jminus(hs=hs)*Jplus(hs=hs)).expand()))
        2 * J_z^(1)
        >>> Jplus(hs=hs).dag() == Jminus(hs=hs)
        True
        >>> Jminus(hs=hs).dag() == Jplus(hs=hs)
        True

    Printers should represent this operator with the default identifier::

        >>> Jz._identifier
        'J_z'

    A custom identifier may be define using `hs`'s `local_identifiers`
    argument.
    """

    _identifier = 'J_z'

    def __init__(self, *, hs):
        super().__init__(hs=hs)

    def _adjoint(self):
        return self

    def _pseudo_inverse(self):
        return PseudoInverse(self) class Jplus(SpinOperator):
    """Raising operator of a spin space

    $\Op{J}_{+} = \Op{J}_x + i \op{J}_y$ is the raising ladder operator
    of a general spin operator acting on a particular :class:`SpinSpace` `hs`
    with well defined spin quantum number $s$.  It's adjoint is the
    lowering operator::

        >>> hs = SpinSpace(1, spin=(1, 2))
        >>> print(ascii(Jplus(hs=hs).adjoint()))
        J_-^(1)

    :class:`Jz`, :class:`Jplus` and :class:`Jminus` satisfy that angular
    momentum commutator algebra, see :class:`Jz`

    Printers should represent this operator with the default identifier::

        >>> Jplus._identifier
        'J_+'

    A custom identifier may be define using `hs`'s `local_identifiers`
    argument.
    """

    _identifier = 'J_+'

    def __init__(self, *, hs):
        super().__init__(hs=hs)

    def _adjoint(self):
        return Jminus(hs=self.space)

    def _pseudo_inverse(self):
        return PseudoInverse(self) class Jminus(SpinOperator):
    """Lowering operator on a spin space

    $\Op{J}_{-} = \Op{J}_x - i \op{J}_y$ is the lowering ladder operator of
    a general spin operator acting on a particular :class:`SpinSpace` `hs`
    with well defined spin quantum number $s$.  It's adjoint is the raising
    operator::

        >>> hs = SpinSpace(1, spin=(1, 2))
        >>> print(ascii(Jminus(hs=hs).adjoint()))
        J_+^(1)

    :class:`Jz`, :class:`Jplus` and :class:`Jminus` satisfy that angular
    momentum commutator algebra, see :class:`Jz`.

    Printers should represent this operator with the default identifier::

        >>> Jminus._identifier
        'J_-'

    A custom identifier may be define using `hs`'s `local_identifiers`
    argument.
    """

    _identifier = 'J_-'

    def __init__(self, *, hs):
        super().__init__(hs=hs)

    def _adjoint(self):
        return Jplus(hs=self.space)

    def _pseudo_inverse(self):
        return PseudoInverse(self) class Destroy(LocalOperator):
    """Bosonic annihilation operator

    It obeys the bosonic commutation relation::

        >>> Destroy(hs=1) * Create(hs=1) - Create(hs=1) * Destroy(hs=1)
        IdentityOperator
        >>> Destroy(hs=1) * Create(hs=2) - Create(hs=2) * Destroy(hs=1)
        ZeroOperator
    """

    _identifier = 'a'
    _dagger = False
    _rx_identifier = re.compile('^[A-Za-z][A-Za-z0-9]*$')

    def __init__(self, *, hs):
        super().__init__(hs=hs)

    def _adjoint(self):
        return Create(hs=self.space)

    def _pseudo_inverse(self):
        return PseudoInverse(self)

    @property
    def identifier(self):
        """The identifier (symbol) that is used when printing the annihilation
        operator. This is identical to the identifier of :class:`Create`. A
        custom identifier for both :class:`Destroy` and :class:`Create` can be
        set through the `local_identifiers` parameter of the associated Hilbert
        space::

            >>> hs_custom = LocalSpace(0, local_identifiers={'Destroy': 'b'})
            >>> Create(hs=hs_custom).identifier
            'b'
            >>> Destroy(hs=hs_custom).identifier
            'b'
        """
        identifier = self._hs._local_identifiers.get(
            self.__class__.__name__, self._hs._local_identifiers.get(
                'Create', self._identifier))
        if not self._rx_identifier.match(identifier):
            raise ValueError(
                "identifier '%s' does not match pattern '%s'"
                % (identifier, self._rx_identifier.pattern))
        return identifier class Create(LocalOperator):
    """Bosonic creation operator

    This is the adjoint of :class:`Destroy`.
    """

    _identifier = 'a'
    _dagger = True
    _rx_identifier = re.compile('^[A-Za-z][A-Za-z0-9]*$')

    def __init__(self, *, hs):
        super().__init__(hs=hs)

    def _adjoint(self):
        return Destroy(hs=self.space)

    def _pseudo_inverse(self):
        return PseudoInverse(self)

    @property
    def identifier(self):
        """The identifier (symbols) that is used when printing the creation
        operator. This is identical to the identifier of :class:`Destroy`"""
        identifier = self._hs._local_identifiers.get(
            self.__class__.__name__, self._hs._local_identifiers.get(
                'Destroy', self._identifier))
        if not self._rx_identifier.match(identifier):
            raise ValueError(
                "identifier '%s' does not match pattern '%s'"
                % (identifier, self._rx_identifier.pattern))
        return identifier class Phase(LocalOperator):
    r"""Unitary "phase" operator

    .. math::

        \op{P}_{\rm hs}(\phi) =
        \exp\left(i \phi \Op{a}_{\rm hs}^\dagger \Op{a}_{\rm hs}\right)

    where :math:`a_{\rm hs}` is the annihilation operator acting on the
    :class:`.LocalSpace` `hs`.

    Args:
        phase (Scalar): the phase $\phi$
        hs (HilbertSpace or int or str): The Hilbert space on which the
            operator acts

    Printers should represent this operator with the default identifier::

        >>> Phase._identifier
        'Phase'

    A custom identifier may be define using `hs`'s `local_identifiers`
    argument.
    """
    _identifier = 'Phase'
    _arg_names = ('phase', )
    _rules = OrderedDict()
    simplifications = [implied_local_space(keys=['hs', ]), match_replace]

    def _adjoint(self):
        return Phase.create(-self.phase.conjugate(), hs=self.space)

    def _pseudo_inverse(self):
        return Phase.create(-self.phase, hs=self.space) class Displace(LocalOperator):
    r"""Unitary coherent displacement operator

    .. math::

        \op{D}_{\rm hs}(\alpha) =
        \exp\left({\alpha \Op{a}_{\rm hs}^\dagger -
                   \alpha^* \Op{a}_{\rm hs}}\right)

    where :math:`\Op{a}_{\rm hs}` is the annihilation operator acting on the
    :class:`.LocalSpace` `hs`.

    Args:
        displacement (Scalar): the displacement amplitude $\alpha$
        hs (HilbertSpace or int or str): The Hilbert space on which the
            operator acts

    Printers should represent this operator with the default identifier::

        >>> Displace._identifier
        'D'

    A custom identifier may be define using `hs`'s `local_identifiers`
    argument.
    """
    _identifier = 'D'
    _nargs = 1
    _arg_names = ('displacement', )
    _rules = OrderedDict()
    simplifications = [implied_local_space(keys=['hs', ]), match_replace]

    def _adjoint(self):
        return Displace.create(-self.displacement, hs=self.space)

    _pseudo_inverse = _adjoint class Squeeze(LocalOperator):
    r"""Unitary squeezing operator

    .. math::

        \Op{S}_{\rm hs}(\eta) =
        \exp {\left( \frac{\eta}{2} {\Op{a}_{\rm hs}^\dagger}^2 -
                     \frac{\eta^*}{2} {\Op{a}_{\rm hs}}^2 \right)}

    where :math:`\Op{a}_{\rm hs}` is the annihilation operator acting on the
    :class:`.LocalSpace` `hs`.

    Args:
        squeezing_factor (Scalar): the squeezing factor $\eta$
        hs (HilbertSpace or int or str): The Hilbert space on which the
            operator acts

    Printers should represent this operator with the default identifier::

        >>> Squeeze._identifier
        'Squeeze'

    A custom identifier may be define using `hs`'s `local_identifiers`
    argument.
    """
    _identifier = "Squeeze"
    _arg_names = ('squeezing_factor', )
    _rules = OrderedDict()
    simplifications = [implied_local_space(keys=['hs', ]), match_replace]

    def _adjoint(self):
        return Squeeze(-self.squeezing_factor, hs=self.space)

    _pseudo_inverse = _adjoint class State(QuantumExpression, metaclass=ABCMeta):
    """Base class for states in a Hilbert space"""

    def _adjoint(self):
        if self.isket:
            return Bra(self)
        else:
            return self.ket

    @property
    def isket(self):
        """Whether the state represents a ket"""
        # We enforce that all sub-classes of Ket except Bra are normal
        # Hilbert space vectors. For example, KetPlus(Bra(...), Bra(...))
        # is disallowed, and is represented as Bra(KetPlus(.., ...)). Thus,
        # all State instances that are not directly an instance of Bra are
        # kets. Name of property 'isket' is taken from QuTiP
        return not isinstance(self, Bra)

    @property
    def isbra(self):
        """Wether the state represents a bra (adjoint ket)"""
        return isinstance(self, Bra)

    @property
    def bra(self):
        """The bra associated with a ket"""
        if self.isbra:
            return self
        else:
            return self._adjoint()

    @property
    def ket(self):
        """The ket associated with a bra"""
        if self.isket:
            return self
        else:
            return self._adjoint()

    def __mul__(self, other):
        # In the State algebra, any adjoint is a Bra
        # instance at the top level. To enforce this, we need some custom
        # __mul__ implementations, in particular for KetIndexedSums

        def check_if_sum(*args):
            return [isinstance(arg, KetIndexedSum) for arg in args]

        if isinstance(other, State):
            isket = (self.isket, other.isket)
            if isket == (True, True):
                # Ket * Ket
                is_sum = check_if_sum(self, other)
                if is_sum == [True, True]:
                    return QuantumIndexedSum.__mul__(self, other)
                elif is_sum == [True, False]:
                    return QuantumIndexedSum.__mul__(self, other)
                elif is_sum == [False, True]:
                    return QuantumIndexedSum.__rmul__(other, self)
                elif is_sum == [False, False]:
                    return TensorKet.create(self, other)
            elif isket == (True, False):
                # Ket * Bra
                is_sum = check_if_sum(self, other.ket)
                if is_sum == [True, True]:
                    other_ket = other.ket.make_disjunct_indices(self)
                    new_term = self.term * Bra(other_ket.term)
                    new_ranges = self.ranges + other_ket.ranges
                    return new_term.__class__._indexed_sum_cls.create(
                        new_term, *new_ranges)
                elif is_sum == [True, False]:
                    return QuantumIndexedSum.__mul__(self, other)
                elif is_sum == [False, True]:
                    new_term = self * Bra(other.ket.term)
                    return new_term.__class__._indexed_sum_cls.create(
                        new_term, *other.ket.ranges)
                elif is_sum == [False, False]:
                    return KetBra.create(self, other.ket)
            elif isket == (False, True):
                # Bra * Ket
                is_sum = check_if_sum(self.ket, other)
                if is_sum == [True, True]:
                    other = other.make_disjunct_indices(self.ket)
                    new_term = Bra(self.ket.term) * other.term
                    new_ranges = self.ket.ranges + other.ranges
                    return new_term.__class__._indexed_sum_cls.create(
                        new_term, *new_ranges)
                elif is_sum == [True, False]:
                    new_term = Bra(self.ket.term) * other
                    return new_term.__class__._indexed_sum_cls.create(
                        new_term, *self.ket.ranges)
                elif is_sum == [False, True]:
                    return QuantumIndexedSum.__rmul__(other, self)
                elif is_sum == [False, False]:
                    return BraKet.create(self.ket, other)
            elif isket == (False, False):
                # Bra * Bra
                return Bra.create(self.ket * other.ket)
        elif isinstance(other, Operator):
            if self.isbra:
                return Bra.create(other.adjoint() * self.ket)
        try:
            return super().__mul__(other)
        except AttributeError:
            return NotImplemented

    def __rmul__(self, other):
        if self.isket and isinstance(other, Operator):
            return OperatorTimesKet.create(other, self)
        elif self.isbra and is_scalar(other):
            return Bra(other.conjugate() * self.ket)
        try:
            return super().__rmul__(other)
        except AttributeError:
            return NotImplemented class BraKet(ScalarExpression, Operation):
    r"""The symbolic inner product between two states

    This mathermatically corresponds to:

    .. math::
        \langle b | k \rangle

    which we define to be linear in the state :math:`k` and anti-linear in
    :math:`b`.

    Args:
        bra (State): The anti-linear state argument. Note that this is *not* a
            :class:`Bra` instance.
        ket (State): The linear state argument.
    """
    _rules = OrderedDict()
    _space = TrivialSpace
    simplifications = [match_replace]

    def __init__(self, bra, ket):
        _check_kets(bra, ket, same_space=True)
        super().__init__(bra, ket)

    @property
    def ket(self):
        """The ket of the braket"""
        return self.operands[1]

    @property
    def bra(self):
        """The bra of the braket (:class:`Bra` instance)"""
        return Bra(self.operands[0])

    def _diff(self, sym):
        bra, ket = self.operands
        return (
            self.__class__.create(bra.diff(sym), ket) +
            self.__class__.create(bra, ket.diff(sym)))

    def _adjoint(self):
        return BraKet.create(*reversed(self.operands))

    def _expand(self):
        b, k = self.operands
        be, ke = b.expand(), k.expand()
        besummands = be.operands if isinstance(be, KetPlus) else (be,)
        kesummands = ke.operands if isinstance(ke, KetPlus) else (ke,)
        return sum(BraKet.create(bes, kes)
                   for bes in besummands for kes in kesummands)

    def _series_expand(self, param, about, order):
        be = self.bra.series_expand(param, about, order)
        ke = self.ket.series_expand(param, about, order)
        return _series_expand_combine_prod(be, ke, order) class KetBra(Operator, Operation):
    """Outer product of two states

    Args:
        ket (State): The left factor in the product
        bra (State): The right factor in the product. Note that this is *not* a
            :class:`Bra` instance.
    """
    _rules = OrderedDict()
    simplifications = [match_replace]

    def __init__(self, ket, bra):
        _check_kets(ket, bra, same_space=True)
        super().__init__(ket, bra)

    @property
    def ket(self):
        """The left factor in the product"""
        return self.operands[0]

    @property
    def bra(self):
        """The co-state right factor in the product

        This is a :class:`Bra` instance (unlike the `bra` given to the
        constructor
        """
        return Bra(self.operands[1])

    @property
    def space(self):
        """The Hilbert space of the states being multiplied"""
        return self.operands[0].space

    def _adjoint(self):
        return KetBra.create(*reversed(self.operands))

    def _pseudo_inverse(self):
        return PseudoInverse(self)

    def _diff(self, sym):
        ket, bra = self.operands
        return (
            self.__class__.create(ket.diff(sym), bra) +
            self.__class__.create(ket, bra.diff(sym)))

    def _expand(self):
        k, b = self.ket, self.bra.ket
        be, ke = b.expand(), k.expand()
        kesummands = ke.operands if isinstance(ke, KetPlus) else (ke,)
        besummands = be.operands if isinstance(be, KetPlus) else (be,)
        res_summands = []
        for (k, b) in cartesian_product(kesummands, besummands):
            res_summands.append(KetBra.create(k, b))
        return OperatorPlus.create(*res_summands)

    def _series_expand(self, param, about, order):
        ke = self.ket.series_expand(param, about, order)
        be = self.bra.series_expand(param, about, order)
        return tuple(ke[k] * be[n - k]
                     for n in range(order + 1) for k in range(n + 1)) class BasisKet(LocalKet, KetSymbol):
    """Local basis state, identified by index or label

    Basis kets are orthornormal, and the :meth:`next` and :meth:`prev` methods
    can be used to move between basis states.

    Args:
        label_or_index: If `str`, the label of the basis state (must be an
            element of `hs.basis_labels`). If `int`, the (zero-based)
            index of the basis state. This works if `hs` has an unknown
            dimension. For a symbolic index, `label_or_index` can be an
            instance of an appropriate subclass of
            :class:`~qnet.algebra.indices.SymbolicLabelBase`
        hs (LocalSpace): The Hilbert space in which the basis is defined

    Raises:
        ValueError: if `label_or_index` is not in the Hilbert space
        TypeError: if `label_or_index` is not of an appropriate type
        .BasisNotSetError: if `label_or_index` is a `str` but no basis is
            defined for `hs`

    Note:
        Basis states that are instantiated via a label or via an index are
        equivalent::

            >>> hs = LocalSpace('tls', basis=('g', 'e'))
            >>> BasisKet('g', hs=hs) == BasisKet(0, hs=hs)
            True
            >>> print(ascii(BasisKet(0, hs=hs)))
            |g>^(tls)

        When instantiating the :class:`BasisKet` via
        :meth:`~qnet.algebra.abstract_algebra.Expression.create`, an integer
        label outside the range of the underlying Hilbert space results in a
        :class:`ZeroKet`::

            >>> BasisKet.create(-1, hs=0)
            ZeroKet
            >>> BasisKet.create(2, hs=LocalSpace('tls', dimension=2))
            ZeroKet
    """

    simplifications = [basis_ket_zero_outside_hs]

    def __init__(self, label_or_index, *, hs):
        hs = ensure_local_space(hs, cls=self._default_hs_cls)
        label, ind = hs._unpack_basis_label_or_index(label_or_index)
        self._index = ind
        super().__init__(label, hs=hs)

    @property
    def args(self):
        """Tuple containing `label_or_index` as its only element."""
        if self.space.has_basis or isinstance(self.label, SymbolicLabelBase):
            return (self.label, )
        else:
            return (self.index, )

    @property
    def index(self):
        """The index of the state in the Hilbert space basis

        >>> hs = LocalSpace('tls', basis=('g', 'e'))
        >>> BasisKet('g', hs=hs).index
        0
        >>> BasisKet('e', hs=hs).index
        1
        >>> BasisKet(1, hs=hs).index
        1

        For a :class:`BasisKet` with an indexed label, this may return a sympy
        expression::

        >>> hs = SpinSpace('s', spin='3/2')
        >>> i = symbols('i', cls=IdxSym)
        >>> lbl = SpinIndex(i/2, hs)
        >>> ket = BasisKet(lbl, hs=hs)
        >>> ket.index
        i/2 + 3/2
        """
        return self._index

    def next(self, n=1):
        """Move up by `n` steps in the Hilbert space::

            >>> hs =  LocalSpace('tls', basis=('g', 'e'))
            >>> ascii(BasisKet('g', hs=hs).next())
            '|e>^(tls)'
            >>> ascii(BasisKet(0, hs=hs).next())
            '|e>^(tls)'

        We can also go multiple steps:

            >>> hs =  LocalSpace('ten', dimension=10)
            >>> ascii(BasisKet(0, hs=hs).next(2))
            '|2>^(ten)'

        An increment that leads out of the Hilbert space returns zero::

            >>> BasisKet(0, hs=hs).next(10)
            ZeroKet

        """
        if isinstance(self.label, SymbolicLabelBase):
            next_label = self.space.next_basis_label_or_index(
                self.label, n)
            return BasisKet(next_label, hs=self.space)
        else:
            try:
                next_index = self.space.next_basis_label_or_index(
                    self.index, n)
                return BasisKet(next_index, hs=self.space)
            except IndexError:
                return ZeroKet

    def prev(self, n=1):
        """Move down by `n` steps in the Hilbert space, cf. :meth:`next`.

        >>> hs =  LocalSpace('3l', basis=('g', 'e', 'r'))
        >>> ascii(BasisKet('r', hs=hs).prev(2))
        '|g>^(3l)'
        >>> BasisKet('r', hs=hs).prev(3)
        ZeroKet
        """
        return self.next(n=-n) class CoherentStateKet(LocalKet):
    """Local coherent state, labeled by a complex amplitude

    Args:
        hs (LocalSpace): The local Hilbert space degree of freedom.
        ampl (Scalar): The coherent displacement amplitude.
    """

    _rx_label = re.compile('^.*$')

    def __init__(self, ampl, *, hs):
        self._ampl = ampl
        super().__init__(ampl, hs=hs)

    @property
    def args(self):
        return (self._ampl, )

    @property
    def ampl(self):
        return self._ampl

    def _diff(self, sym):
        from qnet.algebra.library.fock_operators import Destroy, Create
        hs = self.space
        return (
            (self._ampl * Create(hs=hs) -
             self._ampl.conjugate() * Destroy(hs=hs))
            .diff(sym)
            * self)

    def to_fock_representation(self, index_symbol='n', max_terms=None):
        """Return the coherent state written out as an indexed sum over Fock
        basis states"""
        phase_factor = sympy.exp(
            sympy.Rational(-1, 2) * self.ampl * self.ampl.conjugate())
        if not isinstance(index_symbol, IdxSym):
            index_symbol = IdxSym(index_symbol)
        n = index_symbol
        if max_terms is None:
            index_range = IndexOverFockSpace(n, hs=self._hs)
        else:
            index_range = IndexOverRange(n, 0, max_terms-1)
        term = (
            (self.ampl**n / sympy.sqrt(sympy.factorial(n))) *
            BasisKet(FockIndex(n), hs=self._hs))
        return phase_factor * KetIndexedSum(term, index_range) class KetPlus(State, QuantumPlus):
    """Sum of states"""
    _neutral_element = ZeroKet
    _binary_rules = OrderedDict()
    simplifications = [
        accept_bras, assoc, orderby, collect_summands, match_replace_binary]

    order_key = FullCommutativeHSOrder

    def __init__(self, *operands):
        _check_kets(*operands, same_space=True)
        super().__init__(*operands) class TensorKet(State, QuantumTimes):
    """A tensor product of kets

    Each ket must belong to different degree of freedom (:class:`.LocalSpace`).
    """
    _binary_rules = OrderedDict()
    _neutral_element = TrivialKet
    simplifications = [
        accept_bras, assoc, orderby, filter_neutral, match_replace_binary]

    order_key = FullCommutativeHSOrder

    def __init__(self, *operands):
        _check_kets(*operands, disjunct_space=True)
        super().__init__(*operands)

    @classmethod
    def create(cls, *ops):
        if any(o == ZeroKet for o in ops):
            return ZeroKet
        return super().create(*ops) class ScalarTimesKet(State, ScalarTimesQuantumExpression):
    """Product of a :class:`.Scalar` coefficient and a ket

    Args:
        coeff (Scalar): coefficient
        term (State): the ket that is multiplied
    """
    _rules = OrderedDict()
    simplifications = [match_replace, ]

    @classmethod
    def create(cls, coeff, term):
        if term.isbra:
            scalar_times_ket = coeff.conjugate() * term.ket
            return Bra.create(scalar_times_ket)
        return super().create(coeff, term)

    def __init__(self, coeff, term):
        _check_kets(term)
        super().__init__(coeff, term) class OperatorTimesKet(State, Operation):
    """Product of an operator and a state."""
    _rules = OrderedDict()
    simplifications = [match_replace]

    def __init__(self, operator, ket):
        _check_kets(ket)
        if not operator.space <= ket.space:
            raise SpaceTooLargeError(
                str(operator.space) + " <!= " + str(ket.space))
        super().__init__(operator, ket)

    @property
    def space(self):
        return self.operands[1].space

    @property
    def operator(self):
        return self.operands[0]

    @property
    def ket(self):
        return self.operands[1]

    def _expand(self):
        c, t = self.operands
        ct = c.expand()
        et = t.expand()
        if isinstance(et, KetPlus):
            if isinstance(ct, OperatorPlus):
                return sum((cto * eto
                            for eto in et.operands
                            for cto in ct.operands), ZeroKet)
            else:
                return sum((c * eto for eto in et.operands), ZeroKet)
        elif isinstance(ct, OperatorPlus):
            return sum((cto * et for cto in ct.operands), ZeroKet)
        return ct * et

    def _diff(self, sym):
        return (
            self.operator.diff(sym) * self.ket +
            self.operator * self.ket.diff(sym))

    def _series_expand(self, param, about, order):
        ce = self.operator.series_expand(param, about, order)
        te = self.ket.series_expand(param, about, order)

        return tuple(ce[k] * te[n - k]
                     for n in range(order + 1) for k in range(n + 1)) class TrivialSpace(HilbertSpace, Expression, metaclass=Singleton):
    """The 'nullspace', i.e. a one dimensional Hilbert space, which is a factor
    space of every other Hilbert space.

    This is the Hilbert space of scalars.
    """

    def __hash__(self):
        return hash(self.__class__)

    @property
    def args(self):
        """Empty tuple (no arguments)"""
        return ()

    @property
    def _order_key(self):
        return KeyTuple((-2, '_'))

    @property
    def dimension(self):
        return 1

    @property
    def has_basis(self):
        """True, by definition (the basis is defined as
        :obj:`~qnet.algebra.state_algebra.TrivialKet`)"""
        return True

    @property
    def basis_states(self):
        """Yield the :obj:`~qnet.algebra.state_algebra.TrivialKet`"""
        from qnet.algebra.core.state_algebra import TrivialKet
        yield TrivialKet

    def basis_state(self, index_or_label):
        """Return the basis state with the given index 0 or label "0".

        All other indices or labels raise an exception.

        Raises:
            IndexError: if index is different from 0
            KeyError: if label is differnt from ""
        """
        from qnet.algebra.core.state_algebra import TrivialKet
        if index_or_label in [0, "0"]:
            return TrivialKet
        else:
            if isinstance(index_or_label, int):
                raise IndexError(
                    "No index %d in basis for TrivialSpace" % index_or_label)
            else:
                raise KeyError(
                    "No label %d in basis for TrivialSpace" % index_or_label)

    @property
    def basis_labels(self):
        """The one-element tuple containing the label '0'"""
        return tuple(["0", ])

    def remove(self, other):
        """Removing any Hilbert space from the trivial space yields the trivial
        space again"""
        return self

    def intersect(self, other):
        """The intersection of the trivial space with any other space is only
        the trivial space"""
        return self

    @property
    def local_factors(self):
        """Empty list (the trivial space has no factors)"""
        return ()

    def is_strict_subfactor_of(self, other):
        """The trivial space is a subfactor of any other space (except
        itself)"""
        if other is TrivialSpace:
            return False
        return True

    @property
    def label(self):
        return "null" class SuperOperator(QuantumExpression, metaclass=ABCMeta):
    """Base class for super-operators"""
    def __mul__(self, other):
        if isinstance(other, Operator):
            return SuperOperatorTimesOperator.create(self, other)
        else:
            return super().__mul__(other) class IdentitySuperOperator(SuperOperator, metaclass=Singleton):
    """Neutral element for product of super-operators"""

    _order_index = 2

    @property
    def space(self):
        return TrivialSpace

    @property
    def args(self):
        return tuple()

    def _diff(self, sym):
        return ZeroSuperOperator

    def _adjoint(self):
        return self

    def _expand(self):
        return self class SuperOperatorPlus(QuantumPlus, SuperOperator):
    """A sum of super-operators"""

    _neutral_element = ZeroSuperOperator
    _binary_rules = OrderedDict()
    simplifications = [assoc, orderby, collect_summands, match_replace_binary] class SuperOperatorTimes(QuantumTimes, SuperOperator):
    """Product of super-operators"""

    _neutral_element = IdentitySuperOperator
    _binary_rules = OrderedDict()  # see end of module
    simplifications = [assoc, orderby, filter_neutral, match_replace_binary]

    order_key = SuperCommutativeHSOrder

    @classmethod
    def create(cls, *ops):
        # TODO: Add this functionality to QuantumTimes
        if any(o == ZeroSuperOperator for o in ops):
            return ZeroSuperOperator
        return super().create(*ops) class ScalarTimesSuperOperator(SuperOperator, ScalarTimesQuantumExpression):
    """Product of a :class:`.Scalar` coefficient and a
    :class:`SuperOperator`"""

    def _adjoint(self):
        pass

    _rules = OrderedDict()  # see end of module
    simplifications = [match_replace, ] class SPre(SuperOperator, Operation):
    """Linear pre-multiplication operator

    Acting ``SPre(A)`` on an operator ``B`` just yields the product ``A * B``
    """

    _rules = OrderedDict()  # see end of module
    simplifications = [match_replace, ]

    _order_name = 'A_SPre'  # "SPre" should go before "SPost"

    @property
    def space(self):
        return self.operands[0].space

    def _expand(self):
        oe = self.operands[0].expand()
        if isinstance(oe, OperatorPlus):
            return sum(SPre.create(oet) for oet in oe.operands)
        return SPre.create(oe)

    def _simplify_scalar(self, func):
        return self.create(self.operands[0].simplify_scalar(func=func))

    def _diff(self, sym):
        return SuperOperatorDerivative(self, {sym: 1})

    def _adjoint(self):
        return SPost(self.operands[0]) class SPost(SuperOperator, Operation):
    """Linear post-multiplication operator

        Acting ``SPost(A)`` on an operator ``B`` just yields the reversed
        product ``B * A``.
    """

    _order_index = -1

    _rules = OrderedDict()  # see end of module
    simplifications = [match_replace, ]

    _order_name = 'B_SPost'  # "SPost" should go before "SPre"

    @property
    def space(self):
        return self.operands[0].space

    def _expand(self):
        oe = self.operands[0].expand()
        if isinstance(oe, OperatorPlus):
            return sum(SPost.create(oet) for oet in oe.operands)
        return SPost.create(oe)

    def _simplify_scalar(self, func):
        return self.create(self.operands[0].simplify_scalar(func=func))

    def _diff(self, sym):
        return SuperOperatorDerivative(self, {sym: 1})

    def _adjoint(self):
        return SPre(self.operands[0]) class SuperOperatorTimesOperator(Operator, Operation):
    """Application of a super-operator to an operator

    The result of this operation is(result is an :class:`Operator`
    """

    _rules = OrderedDict()  # see end of module
    simplifications = [match_replace, ]

    def __init__(self, sop, op):
        assert isinstance(sop, SuperOperator)
        assert isinstance(op, Operator)
        super().__init__(sop, op)

    @property
    def space(self):
        return self.sop.space * self.op.space

    @property
    def sop(self):
        return self.operands[0]

    @property
    def op(self):
        return self.operands[1]

    def _expand(self):
        sop, op = self.operands
        sope, ope = sop.expand(), op.expand()
        if isinstance(sope, SuperOperatorPlus):
            sopet = sope.operands
        else:
            sopet = (sope, )
        if isinstance(ope, OperatorPlus):
            opet = ope.operands
        else:
            opet = (ope, )

        return sum(st * ot for st in sopet for ot in opet)

    def _series_expand(self, param, about, order):
        sop, op = self.sop, self.op
        ope = op.series_expand(param, about, order)
        return tuple(sop * opet for opet in ope)

    def _simplify_scalar(self, func):
        sop, op = self.sop, self.op
        return sop.simplify_scalar(func=func) * op.simplify_scalar(func=func)

    def _adjoint(self):
        return Adjoint(self)

    def _pseudo_inverse(self):
        return PseudoInverse(self)

    def _diff(self, sym):
        return self.sop.diff(sym) * self.op + self.sop * self.op.diff(sym) class ZeroSuperOperator(SuperOperator, metaclass=Singleton):
    """Neutral element for sum of super-operators"""

    _order_index = 2

    @property
    def space(self):
        return TrivialSpace

    @property
    def args(self):
        return tuple()

    def _diff(self, sym):
        return self

    def _adjoint(self):
        return self

    def _expand(self):
        return self
import re
    import qutip
from functools import reduce
from sympy import symbols
from sympy.utilities.lambdify import lambdify
from scipy.sparse import csr_matrix
from numpy import (
    diag as np_diag, arange, cos as np_cos, sin as np_sin)
from ..algebra.core.scalar_algebra import ScalarValue, is_scalar
from ..algebra.core.exceptions import AlgebraError
from ..algebra.core.circuit_algebra import SLH, move_drive_to_H
from ..algebra.core.abstract_algebra import Operation
from ..algebra.core.operator_algebra import (
    Operator, IdentityOperator, ZeroOperator, LocalOperator, LocalSigma,
    OperatorPlus, OperatorTimes, ScalarTimesOperator,
    Adjoint, PseudoInverse, OperatorTrace, NullSpaceProjector)
from ..algebra.library.spin_algebra import Jz, Jplus, Jminus
from ..algebra.library.fock_operators import (
    Destroy, Create, Phase, Displace, Squeeze)
from ..algebra.core.state_algebra import (
    State, BraKet, KetBra, BasisKet, CoherentStateKet, KetPlus, TensorKet,
    ScalarTimesKet, OperatorTimesKet)
from ..algebra.core.hilbert_space_algebra import TrivialSpace
from ..algebra.core.super_operator_algebra import (
    SuperOperator, IdentitySuperOperator, SuperOperatorPlus,
    SuperOperatorTimes, ScalarTimesSuperOperator, SPre, SPost,
    SuperOperatorTimesOperator, ZeroSuperOperator)
            from scipy.linalg import pinv
            from scipy.linalg import svd
        ``H1`` are of type `qutip.Qobj`, and ``eps_t`` is either a string
        (``convert_as='str'``) or a function (``convert_as='pyfunc'``)

    Raises:
        AlgebraError: If the Hilbert space (`slh.space` or `full_space`) is
            invalid for numerical conversion
    """
    if full_space:
        if not full_space >= slh.space:
            raise AlgebraError("full_space="+str(full_space)+" needs to "
                               "at least include slh.space = "+str(slh.space))
    else:
        full_space = slh.space
    if full_space == TrivialSpace:
        raise AlgebraError(
            "Cannot convert SLH object in TrivialSpace. "
            "You may pass a non-trivial `full_space`")
    slh = move_drive_to_H(slh)
    if time_symbol is None:
        H = convert_to_qutip(slh.H, full_space=full_space)
        Ls = []
        for L in slh.Ls:
            if is_scalar(L):
                L = L * IdentityOperator
            L_qutip = convert_to_qutip(L, full_space=full_space)
            if L_qutip.norm('max') > 0:
                Ls.append(L_qutip)
    else:
        H = _time_dependent_to_qutip(slh.H, full_space, time_symbol,
                                     convert_as)
        Ls = []
        for L in slh.Ls:
            if is_scalar(L):
                L = L * IdentityOperator
            L_qutip = _time_dependent_to_qutip(L, full_space, time_symbol,
                                               convert_as)
            Ls.append(L_qutip)
    return H, Ls


def _convert_local_operator_to_qutip(expr, full_space, mapping):
    """Convert a LocalOperator instance to qutip"""
    n = full_space.dimension
    if full_space != expr.space:
        all_spaces = full_space.local_factors
        own_space_index = all_spaces.index(expr.space)
        return qutip.tensor(
            *([qutip.qeye(s.dimension)
               for s in all_spaces[:own_space_index]] +
              [convert_to_qutip(expr, expr.space, mapping=mapping), ] +
              [qutip.qeye(s.dimension)
               for s in all_spaces[own_space_index + 1:]])
        )
    if isinstance(expr, Create):
        return qutip.create(n)
    elif isinstance(expr, Jz):
        return qutip.jmat((expr.space.dimension-1)/2., "z")
    elif isinstance(expr, Jplus):
        return qutip.jmat((expr.space.dimension-1)/2., "+")
    elif isinstance(expr, Jminus):
        return qutip.jmat((expr.space.dimension-1)/2., "-")
    elif isinstance(expr, Destroy):
        return qutip.destroy(n)
    elif isinstance(expr, Phase):
        arg = complex(expr.operands[1]) * arange(n)
        d = np_cos(arg) + 1j * np_sin(arg)
        return qutip.Qobj(np_diag(d))
    elif isinstance(expr, Displace):
        alpha = expr.operands[1]
        return qutip.displace(n, alpha)
    elif isinstance(expr, Squeeze):
        eta = expr.operands[1]
        return qutip.displace(n, eta)
    elif isinstance(expr, LocalSigma):
        j = expr.j
        k = expr.k
        if isinstance(j, str):
            j = expr.space.basis_labels.index(j)
        if isinstance(k, str):
            k = expr.space.basis_labels.index(k)
        ket = qutip.basis(n, j)
        bra = qutip.basis(n, k).dag()
        return ket * bra
    else:
        raise ValueError("Cannot convert '%s' of type %s"
                         % (str(expr), type(expr)))


def _convert_operator_operation_to_qutip(expr, full_space, mapping):
    if isinstance(expr, OperatorPlus):
        return sum((convert_to_qutip(op, full_space, mapping=mapping)
                    for op in expr.operands), 0)
    elif isinstance(expr, OperatorTimes):
        # if any factor acts non-locally, we need to expand distributively.
        if any(len(op.space) > 1 for op in expr.operands):
            se = expr.expand()
            if se == expr:
                raise ValueError("Cannot represent as QuTiP object: {!s}"
                                 .format(expr))
            return convert_to_qutip(se, full_space, mapping=mapping)
        all_spaces = full_space.local_factors
        by_space = []
        ck = 0
        for ls in all_spaces:
            # group factors by associated local space
            ls_ops = [convert_to_qutip(o, o.space, mapping=mapping)
                      for o in expr.operands if o.space == ls]
            if len(ls_ops):
                # compute factor associated with local space
                by_space.append(reduce(lambda a, b: a * b, ls_ops, 1))
                ck += len(ls_ops)
            else:
                # if trivial action, take identity matrix
                by_space.append(qutip.qeye(ls.dimension))
        assert ck == len(expr.operands)
        # combine local factors in tensor product
        return qutip.tensor(*by_space)
    elif isinstance(expr, Adjoint):
        return convert_to_qutip(qutip.dag(expr.operands[0]), full_space,
                                mapping=mapping)

    elif isinstance(expr, ScalarTimesOperator):