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 ScalarTimesQuantumExpression(
        QuantumExpression, Operation, metaclass=ABCMeta):
    """Product of a :class:`.Scalar` and a :class:`QuantumExpression`"""

    @classmethod
    def create(cls, coeff, term):
        from qnet.algebra.core.scalar_algebra import Scalar, ScalarValue
        if not isinstance(coeff, Scalar):
            coeff = ScalarValue.create(coeff)
        return super().create(coeff, term)

    def __init__(self, coeff, term):
        from qnet.algebra.core.scalar_algebra import Scalar, ScalarValue
        if not isinstance(coeff, Scalar):
            coeff = ScalarValue.create(coeff)
        self._order_coeff = coeff
        self._order_args = KeyTuple([term._order_key])
        super().__init__(coeff, term)

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

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

    def _substitute(self, var_map, safe=False):
        st = self.term.substitute(var_map)
        if isinstance(self.coeff, sympy.Basic):
            svar_map = {k: v for k, v in var_map.items()
                        if not isinstance(k, Expression)}
            sc = self.coeff.subs(svar_map)
        else:
            sc = substitute(self.coeff, var_map)
        if safe:
            return self.__class__(sc, st)
        else:
            return sc * st

    @property
    def free_symbols(self):
        return self.coeff.free_symbols | self.term.free_symbols

    def _adjoint(self):
        return self.coeff.conjugate() * self.term.adjoint()

    @property
    def _order_key(self):
        from qnet.printing.asciiprinter import QnetAsciiDefaultPrinter
        ascii = QnetAsciiDefaultPrinter().doprint
        t = self.term._order_key
        try:
            c = abs(float(self.coeff))  # smallest coefficients first
        except (ValueError, TypeError):
            c = float('inf')
        return KeyTuple(t[:2] + (c,) + t[3:] + (ascii(self.coeff),))

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

    def _expand(self):
        return _expand_product(self.operands)

    def _series_expand(self, param, about, order):
        ce = self.coeff.series_expand(param, about, order)
        te = self.term.series_expand(param, about, order)
        return _series_expand_combine_prod(ce, te, order)

    def _diff(self, sym):
        c, t = self.operands
        return c.diff(sym) * t + c * t.diff(sym)

    def _simplify_scalar(self, func):
        coeff, term = self.operands
        try:
            if isinstance(coeff.val, sympy.Basic):
                coeff = func(coeff)
        except AttributeError:
            # coeff is not a SymPy ScalarValue; leave it unchanged
            pass
        return coeff * term.simplify_scalar(func=func)

    def __complex__(self):
        if self.term is self.__class__._one:
            return complex(self.coeff)
        return NotImplemented

    def __float__(self):
        if self.term is self.__class__._one:
            return float(self.coeff)
        return NotImplemented class QuantumExpression(Expression, metaclass=ABCMeta):
    """Base class for expressions associated with a Hilbert space"""

    _zero = None  # The neutral element for addition
    _one = None  # The neutral element for multiplication
    _base_cls = None  # The most general class we can add / multiply
    _scalar_times_expr_cls = None   # class for multiplication with scalar
    _plus_cls = None  # class for internal addition
    _times_cls = None  # class for internal multiplication
    _adjoint_cls = None  # class for the adjoint
    _indexed_sum_cls = None  # class for indexed sum
    _default_hs_cls = None  # class for implicit Hilbert spaces (str, int)
    # _default_hs_cls is set by `init_algebra`

    _order_index = 0  # index of "order group": things that should go together
    _order_coeff = 1  # scalar prefactor
    _order_name = None

    def __init__(self, *args, **kwargs):
        self._order_args = KeyTuple([
            arg._order_key if hasattr(arg, '_order_key') else arg
            for arg in args])
        self._order_kwargs = KeyTuple([
            KeyTuple([
                key, val._order_key if hasattr(val, '_order_key') else val])
            for (key, val) in sorted(kwargs.items())])
        super().__init__(*args, **kwargs)

    @property
    def is_zero(self):
        """Check whether the expression is equal to zero.

        Specifically, this checks whether the expression is equal to the
        neutral element for the addition within the algebra. This does not
        generally imply equality with a scalar zero:

        >>> ZeroOperator.is_zero
        True
        >>> ZeroOperator == 0
        False
        """
        return self == self._zero

    @property
    def _order_key(self):
        return KeyTuple([
            self._order_index, self._order_name or self.__class__.__name__,
            self._order_coeff, self._order_args, self._order_kwargs])

    @property
    @abstractmethod
    def space(self):
        """The :class:`.HilbertSpace` on which the operator acts
        non-trivially"""
        raise NotImplementedError(self.__class__.__name__)

    def adjoint(self):
        """The Hermitian adjoint of the Expression"""
        return self._adjoint()

    def dag(self):
        """Alias for :meth:`adjoint`"""
        return self._adjoint()

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

    def expand(self):
        """Expand out distributively all products of sums.

        Note:
            This does not expand out sums of scalar coefficients. You may use
            :meth:`simplify_scalar` for this purpose.
        """
        return self._expand()

    def _expand(self):
        return self

    def simplify_scalar(self, func=sympy.simplify):
        """Simplify all scalar symbolic (SymPy) coefficients by appyling `func`
        to them"""
        return self._simplify_scalar(func=func)

    def _simplify_scalar(self, func):
        return self

    def diff(self, sym: Symbol, n: int = 1, expand_simplify: bool = True):
        """Differentiate by scalar parameter `sym`.

        Args:
            sym: What to differentiate by.
            n: How often to differentiate
            expand_simplify: Whether to simplify the result.

        Returns:
            The n-th derivative.
        """
        if not isinstance(sym, sympy.Basic):
            raise TypeError("%s needs to be a Sympy symbol" % sym)
        if sym.free_symbols.issubset(self.free_symbols):
            # QuantumDerivative.create delegates internally to _diff (the
            # explicit non-trivial derivative). Using `create` gives us free
            # caching
            deriv = QuantumDerivative.create(self, derivs={sym: n}, vals=None)
            if not deriv.is_zero and expand_simplify:
                deriv = deriv.expand().simplify_scalar()
            return deriv
        else:
            # the "issubset" of free symbols is a sufficient, but not a
            # necessary condition; if `sym` is non-atomic, determining whether
            # `self` depends on `sym` is not completely trivial (you'd have to
            # substitute with a Dummy)
            return self.__class__._zero

    @abstractmethod
    def _diff(self, sym):
        raise NotImplementedError()

    def series_expand(
            self, param: Symbol, about, order: int) -> tuple:
        r"""Expand the expression as a truncated power series in a
        scalar parameter.

        When expanding an expr for a parameter $x$ about the point $x_0$ up to
        order $N$, the resulting coefficients $(c_1, \dots, c_N)$ fulfill

        .. math::

            \text{expr} = \sum_{n=0}^{N} c_n (x - x_0)^n + O(N+1)

        Args:
            param: Expansion parameter $x$
            about (Scalar): Point $x_0$ about which to expand
            order: Maximum order $N$ of expansion (>= 0)

        Returns:
            tuple of length ``order + 1``, where the entries are the
            expansion coefficients, $(c_0, \dots, c_N)$.

        Note:
            The expansion coefficients are
            "type-stable", in that they share a common base class with the
            original expression. In particular, this applies to "zero"
            coefficients::

                >>> expr = KetSymbol("Psi", hs=0)
                >>> t = sympy.symbols("t")
                >>> assert expr.series_expand(t, 0, 1) == (expr, ZeroKet)
        """
        expansion = self._series_expand(param, about, order)
        # _series_expand is generally not "type-stable", so we continue to
        # ensure the type-stability
        res = []
        for v in expansion:
            if v == 0 or v.is_zero:
                v = self._zero
            elif v == 1:
                v = self._one
            assert isinstance(v, self._base_cls)
            res.append(v)
        return tuple(res)

    def _series_expand(self, param, about, order):
        expr = self
        result = [_evaluate_at(expr, param, about)]
        for n in range(1, order+1):
            if not expr.is_zero:
                expr = expr.diff(param)
                result.append(
                    _evaluate_at(expr, param, about) / sympy.factorial(n))
            else:
                result.append(expr)
        return tuple(result)

    def __add__(self, other):
        if not isinstance(other, self._base_cls):
            try:
                other = self.__class__._one * other
            except TypeError:
                pass
        if isinstance(other, self.__class__._base_cls):
            return self.__class__._plus_cls.create(self, other)
        else:
            return NotImplemented

    def __radd__(self, other):
        # addition is assumed to be commutative
        return self.__add__(other)

    def __mul__(self, other):
        from qnet.algebra.core.scalar_algebra import is_scalar, ScalarValue
        if not isinstance(other, self._base_cls):
            if is_scalar(other):
                other = ScalarValue.create(other)
                # if other was an ScalarExpression, the conversion above leaves
                # it unchanged
                return self.__class__._scalar_times_expr_cls.create(
                    other, self)
        if isinstance(other, self.__class__._base_cls):
            return self.__class__._times_cls.create(self, other)
        else:
            return NotImplemented

    def __rmul__(self, other):
        # multiplication with scalar is assumed to be commutative, but any
        # other multiplication is not
        from qnet.algebra.core.scalar_algebra import is_scalar
        if is_scalar(other):
            return self.__mul__(other)
        else:
            return NotImplemented

    def __sub__(self, other):
        return self + (-1) * other

    def __rsub__(self, other):
        return (-1) * self + other

    def __neg__(self):
        return (-1) * self

    def __truediv__(self, other):
        try:
            factor = _sympyOne / other
            return self * factor
        except TypeError:
            try:
                return super().__rmul__(other)
            except AttributeError:
                return NotImplemented

    def __pow__(self, other):
        if other == 0:
            return self._one
        elif other == 1:
            return self
        else:
            try:
                other_is_int = (other == int(other))
            except TypeError:
                other_is_int = False
            if other_is_int:
                if other > 1:
                    return self.__class__._times_cls.create(
                        *[self for _ in range(other)])
                elif other < 1:
                    return 1 / self**(-other)
                else:
                    raise ValueError("Invalid exponent %r" % other)
            else:
                return NotImplemented class QuantumSymbol(QuantumExpression, metaclass=ABCMeta):
    """Symbolic element of an algebra

    Args:
        label (str or SymbolicLabelBase): Label for the symbol
        sym_args (Scalar): optional scalar arguments. With zero `sym_args`, the
            resulting symbol is a constant. With one or more `sym_args`, it
            becomes a function.
        hs (HilbertSpace, str, int, tuple): the Hilbert space associated with
            the symbol. If a `str` or an `int`, an implicit (sub-)instance of
            :class:`LocalSpace` with a corresponding label will be created, or,
            for a tuple of `str` or `int`, a :class:`ProducSpace. The type of
            the implicit Hilbert space is set by :func:`.init_algebra`.
    """
    _rx_label = re.compile('^[A-Za-z][A-Za-z0-9]*(_[A-Za-z0-9().+-=]+)?$')

    def __init__(self, label, *sym_args, hs):
        from qnet.algebra.core.scalar_algebra import ScalarValue
        self._label = label
        sym_args = [ScalarValue.create(arg) for arg in sym_args]
        self._sym_args = tuple(sym_args)
        if isinstance(label, str):
            if not self._rx_label.match(label):
                raise ValueError(
                    "label '%s' does not match pattern '%s'"
                    % (label, self._rx_label.pattern))
        elif isinstance(label, SymbolicLabelBase):
            self._label = label
        else:
            raise TypeError(
                "type of label must be str or SymbolicLabelBase, not %s"
                % type(label))
        if isinstance(hs, (str, int)):
            hs = self._default_hs_cls(hs)
        elif isinstance(hs, tuple):
            hs = ProductSpace.create(*[self._default_hs_cls(h) for h in hs])
        self._hs = hs
        super().__init__(label, *sym_args, hs=hs)

    @property
    def label(self):
        """Label of the symbol"""
        return self._label

    @property
    def args(self):
        """Tuple of positional arguments, consisting of the label and possible
        `sym_args`"""
        return (self.label, ) + self._sym_args

    @property
    def kwargs(self):
        """Dict of keyword arguments, containing only `hs`"""
        return {'hs': self._hs}

    @property
    def sym_args(self):
        """Tuple of scalar arguments of the symbol"""
        return self._sym_args

    @property
    def space(self):
        return self._hs

    def _diff(self, sym):
        if all([arg.diff(sym).is_zero for arg in self.sym_args]):
            # This includes the case where sym_args is empty
            return self.__class__._zero
        else:
            return self.__class__._derivative_cls(self, derivs={sym: 1})

    def _series_expand(self, param, about, order):
        if len(self._sym_args) == 0:
            return (self, ) + (0, ) * order
        else:
            # QuantumExpression._series_expand will return the abstract Taylor
            # series
            return super()._series_expand(param, about, order)

    def _simplify_scalar(self, func):
        simplified_sym_args = [func(sym) for sym in self._sym_args]
        return self.create(self.label, *simplified_sym_args, hs=self.space)

    def _expand(self):
        return self

    @property
    def free_symbols(self):
        try:
            res = self.label.free_symbols
            # TODO: anywhere else there are symbolic labels, symbols from the
            # labels should be included in the free_symbols, too
        except AttributeError:
            res = set()
        res.update(self._hs.free_symbols)
        return res.union(*[sym.free_symbols for sym in self.sym_args])

    def _adjoint(self):
        return self.__class__._adjoint_cls(self) class QuantumOperation(QuantumExpression, Operation, metaclass=ABCMeta):
    """Base class for operations on quantum expression

    These are operations on quantum expressions within the same fundamental
    set."""

    # "same fundamental set" means all operandas are instances of _base_cls
    # Operations that involve objects from different sets (e.g.,
    # OperatorTimesKet) should directly subclass from QuantumExpression and
    # Operation

    _order_index = 1  # Operations are printed after "atomic" Expressions

    def __init__(self, *operands, **kwargs):
        for o in operands:
            assert isinstance(o, self.__class__._base_cls)
        op_spaces = [o.space for o in operands]
        self._space = ProductSpace.create(*op_spaces)
        super().__init__(*operands, **kwargs)

    @property
    def space(self):
        """Hilbert space of the operation result"""
        return self._space

    def _simplify_scalar(self, func):
        simplified_operands = []
        operands_have_changed = False
        for op in self.operands:
            new_op = op.simplify_scalar(func=func)
            simplified_operands.append(new_op)
            if new_op is not op:
                operands_have_changed = True
        if operands_have_changed:
            return self.create(*simplified_operands, **self.kwargs)
        else:
            return self class QuantumPlus(QuantumOperation, metaclass=ABCMeta):
    """General implementation of addition of quantum expressions"""
    order_key = FullCommutativeHSOrder
    _neutral_element = None

    def __init__(self, *operands, **kwargs):
        if len(operands) <= 1:
            raise TypeError(
                "%s requires at least two operands" % self.__class__.__name__)
        super().__init__(*operands, **kwargs)

    def _expand(self):
        summands = [o.expand() for o in self.operands]
        return self.__class__._plus_cls.create(*summands)

    def _series_expand(self, param, about, order):
        tuples = (o.series_expand(param, about, order) for o in self.operands)
        res = (self.__class__._plus_cls.create(*tels) for tels in zip(*tuples))
        return res

    def _diff(self, sym):
        return sum([o.diff(sym) for o in self.operands], self.__class__._zero)

    def _adjoint(self):
        return self.__class__._plus_cls(*[o.adjoint() for o in self.operands]) class QuantumTimes(QuantumOperation, metaclass=ABCMeta):
    """General implementation of product of quantum expressions"""
    order_key = DisjunctCommutativeHSOrder
    _neutral_element = None

    def __init__(self, *operands, **kwargs):
        if len(operands) <= 1:
            raise TypeError(
                "%s requires at least two operands" % self.__class__.__name__)
        super().__init__(*operands, **kwargs)

    def factor_for_space(self, spc):
        """Return a tuple of two products, where the first product contains the
        given Hilbert space, and the second product is disjunct from it."""
        if spc == TrivialSpace:
            ops_on_spc = [
                o for o in self.operands if o.space is TrivialSpace]
            ops_not_on_spc = [
                o for o in self.operands if o.space > TrivialSpace]
        else:
            ops_on_spc = [
                o for o in self.operands if (o.space & spc) > TrivialSpace]
            ops_not_on_spc = [
                o for o in self.operands if (o.space & spc) is TrivialSpace]
        return (
            self.__class__._times_cls.create(*ops_on_spc),
            self.__class__._times_cls.create(*ops_not_on_spc))

    def _expand(self):
        return _expand_product(self.operands)

    def _series_expand(self, param, about, order):
        assert len(self.operands) > 1
        cfirst = self.operands[0].series_expand(param, about, order)
        if len(self.operands) > 2:
            crest = (
                self.__class__(*self.operands[1:])
                .series_expand(param, about, order))
        else:
            # a single remaining operand needs to be treated explicitly because
            # we didn't use `create` for the `crest` above, for efficiency
            crest = self.operands[1].series_expand(param, about, order)
        return _series_expand_combine_prod(cfirst, crest, order)

    def _diff(self, sym):
        assert len(self.operands) > 1
        first = self.operands[0]
        rest = self.__class__._times_cls.create(*self.operands[1:])
        return first.diff(sym) * rest + first * rest.diff(sym)

    def _adjoint(self):
        return self.__class__._times_cls.create(
                *[o.adjoint() for o in reversed(self.operands)]) class QuantumAdjoint(SingleQuantumOperation, metaclass=ABCMeta):
    """Base class for adjoints of quantum expressions"""

    def _expand(self):
        eo = self.operand.expand()
        if isinstance(eo, self.__class__._plus_cls):
            summands = [eoo.adjoin() for eoo in eo.operands]
            return self.__class__._plus_cls.create(*summands)
        return eo.adjoint()

    def _diff(self, sym):
        return self.__class__.create(self.operands[0].diff(sym))

    def _adjoint(self):
        return self.operand class QuantumDerivative(SingleQuantumOperation):
    r"""Symbolic partial derivative

    .. math::

        \frac{\partial^n}{\partial x_1^{n_1} \dots \partial x_N^{n_N}}
            A(x_1, \dots, x_N); \qquad
        \text{with} \quad n = \sum_i n_i

    Alternatively, if `vals` is given, a symbolic representation of the
    derivative (partially) evaluated at a specific point.

    .. math::

        \left.\frac{\partial^n}{\partial x_1^{n_1} \dots \partial x_N^{n_N}}
            A(x_1, \dots, x_N) \right\vert_{x_1=v_1, \dots}

    Args:
        op (QuantumExpression): the expression $A(x_1, \dots, x_N)$ that is
            being derived
        derivs (dict): a map of symbols $x_i$ to the order $n_i$ of the
            derivate with respect to that symbol
        vals (dict or None): If not ``None``,  a map of symbols $x_i$ to values
            $v_i$ for the point at which the derivative should be evaluated.

    Note:
        :class:`QuantumDerivative` is intended to be instantiated only inside
        the :meth:`_diff` method of a :class:`QuantumExpression`, for
        expressions that depend on scalar arguments in an unspecified way.
        Generally, if a derivative can be calculated explicitly, the explicit
        form is preferred over the abstract :class:`QuantumDerivative`.
    """
    simplifications = [derivative_via_diff, ]  # create -> ._diff
    # *Any* invocations of `create` will directly return the result of
    # `derivative_via_diff` (but with caching)

    @classmethod
    def create(cls, op, *, derivs, vals=None):
        """Instantiate the derivative by repeatedly calling
        the :meth:`~QuantumExpression._diff` method of `op` and evaluating the
        result at the given `vals`.
        """
        # To ensure stable ordering in Expression._get_instance_key, we explicitly
        # convert `derivs` and `vals` to a tuple structure with a custom sorting key.
        if not isinstance(derivs, tuple):
            derivs = cls._dict_to_ordered_tuple(dict(derivs))
        if not (isinstance(vals, tuple) or vals is None):
            vals = cls._dict_to_ordered_tuple(dict(vals))
        return super().create(op, derivs=derivs, vals=vals)

    def __init__(self, op, *, derivs, vals=None):
        self._derivs = defaultdict(int)
        self._n = 0
        syms = []
        for sym, n in dict(derivs).items():
            if not isinstance(sym, sympy.Basic):
                raise TypeError("%s needs to be a Sympy symbol" % sym)
            syms.append(sym)
            n = int(n)
            if n <= 0:
                raise ValueError(
                    "Derivative wrt symbol %s must be be for an order "
                    "greater than zero, not %s" % (sym, n))
            assert n > 0
            self._n += n
            self._derivs[sym] += n
        self._syms = set(syms)
        self._vals = {}
        if vals is not None:
            for sym, val in dict(vals).items():
                if sym not in self._syms:
                    raise ValueError(
                        "Derivative can only be evaluated for a symbol if "
                        "the derivative is with respect to that symbol; "
                        "not %s" % sym)
                self._vals[sym] = val

        self._ordered_derivs = self._dict_to_ordered_tuple(self._derivs)
        self._ordered_vals = self._dict_to_ordered_tuple(self._vals)
        # Expression._get_instance_key wouldn't work with mutable dicts
        super().__init__(
            op, derivs=self._ordered_derivs, vals=self._ordered_vals)

    @property
    def kwargs(self):
        """Keyword arguments for the instantiation of the derivative"""
        return OrderedDict([
            ('derivs', self._ordered_derivs),
            ('vals', self._ordered_vals or None)])

    @property
    def minimal_kwargs(self):
        """Minimal keyword arguments for the instantiation of the derivative
        (excluding defaults)"""
        if len(self._vals) == 0:
            return {'derivs': self._ordered_derivs}
        else:
            return self.kwargs

    @staticmethod
    def _dict_to_ordered_tuple(d):
        from qnet.printing.asciiprinter import QnetAsciiDefaultPrinter
        sort_key = QnetAsciiDefaultPrinter().doprint  # arbitrary, but stable
        return tuple([(s, d[s]) for s in sorted(d.keys(), key=sort_key)])

    def evaluate_at(self, vals):
        """Evaluate the derivative at a specific point"""
        new_vals = self._vals.copy()
        new_vals.update(vals)
        return self.__class__(self.operand, derivs=self._derivs, vals=new_vals)

    @property
    def derivs(self):
        """Mapping of symbols to the order of the derivative with respect to
        that symbol. Keys are ordered alphanumerically."""
        return OrderedDict(self._ordered_derivs)

    @property
    def syms(self):
        """Set of symbols with respect to which the derivative is taken"""
        return set(self._syms)

    @property
    def vals(self):
        """Mapping of symbols to values for which the derivative is to be
        evaluated. Keys are ordered alphanumerically."""
        return OrderedDict(self._ordered_vals)

    @property
    def free_symbols(self):
        """Set of free SymPy symbols contained within the expression."""
        if self._free_symbols is None:
            if len(self._vals) == 0:
                self._free_symbols = self.operand.free_symbols
            else:
                dummy_map = {}
                for sym in self._vals.keys():
                    dummy_map[sym] = sympy.Dummy()
                # bound symbols may not be atomic, so we have to replace them
                # with dummies
                self._free_symbols = {
                    sym for sym
                    in self.operand.substitute(dummy_map).free_symbols
                    if not isinstance(sym, sympy.Dummy)}
                for val in self._vals.values():
                    self._free_symbols.update(val.free_symbols)
        return self._free_symbols

    @property
    def bound_symbols(self):
        """Set of Sympy symbols that are eliminated by evaluation."""
        if self._bound_symbols is None:
            res = set()
            self._bound_symbols = res.union(
                *[sym.free_symbols for sym in self._vals.keys()])
        return self._bound_symbols

    @property
    def n(self):
        """The total order of the derivative.

        This is the sum of the order values in :attr:`derivs`
        """
        return self._n

    def _diff(self, sym):
        if sym in self._vals.keys():
            return self.__class__._zero
        else:
            if not isinstance(sym, sympy.Basic):
                raise TypeError("%s must be a Sympy symbol" % sym)
            if sym in self._vals.values():
                return self.__class__(self, derivs={sym: 1})
            else:
                derivs = self._derivs.copy()
                derivs[sym] += 1
                return self.__class__(
                    self.operand, derivs=derivs, vals=self._vals)

    def _adjoint(self):
        return self.__class__(
            self.operand.adjoint(), derivs=self._derivs, vals=self._vals) def assoc(cls, ops, kwargs):
    """Associatively expand out nested arguments of the flat class.
    E.g.::

        >>> class Plus(Operation):
        ...     simplifications = [assoc, ]
        >>> Plus.create(1,Plus(2,3))
        Plus(1, 2, 3)
    """
    expanded = [(o,) if not isinstance(o, cls) else o.operands for o in ops]
    return sum(expanded, ()), kwargs def filter_neutral(cls, ops, kwargs):
    """Remove occurrences of a neutral element from the argument/operand list,
    if that list has at least two elements.  To use this, one must also specify
    a neutral element, which can be anything that allows for an equality check
    with each argument.  E.g.::

        >>> class X(Operation):
        ...     _neutral_element = 1
        ...     simplifications = [filter_neutral, ]
        >>> X.create(2,1,3,1)
        X(2, 3)
    """
    c_n = cls._neutral_element
    if len(ops) == 0:
        return c_n
    fops = [op for op in ops if c_n != op]  # op != c_n does NOT work
    if len(fops) > 1:
        return fops, kwargs
    elif len(fops) == 1:
        # the remaining operand is the single non-trivial one
        return fops[0]
    else:
        # the original list of operands consists only of neutral elements
        return ops[0] def match_replace(cls, ops, kwargs):
    """Match and replace a full operand specification to a function that
    provides a replacement for the whole expression
    or raises a :exc:`.CannotSimplify` exception.
    E.g.

    First define an operation::

        >>> class Invert(Operation):
        ...     _rules = OrderedDict()
        ...     simplifications = [match_replace, ]

    Then some _rules::

        >>> A = wc("A")
        >>> A_float = wc("A", head=float)
        >>> Invert_A = pattern(Invert, A)
        >>> Invert._rules.update([
        ...     ('r1', (pattern_head(Invert_A), lambda A: A)),
        ...     ('r2', (pattern_head(A_float), lambda A: 1./A)),
        ... ])

    Check rule application::

        >>> print(srepr(Invert.create("hallo")))  # matches no rule
        Invert('hallo')
        >>> Invert.create(Invert("hallo"))        # matches first rule
        'hallo'
        >>> Invert.create(.2)                     # matches second rule
        5.0

    A pattern can also have the same wildcard appear twice::

        >>> class X(Operation):
        ...     _rules = {
        ...         'r1': (pattern_head(A, A), lambda A: A),
        ...     }
        ...     simplifications = [match_replace, ]
        >>> X.create(1,2)
        X(1, 2)
        >>> X.create(1,1)
        1

    """
    expr = ProtoExpr(ops, kwargs)
    if LOG:
        logger = logging.getLogger('QNET.create')
    for key, rule in cls._rules.items():
        pat, replacement = rule
        match_dict = match_pattern(pat, expr)
        if match_dict:
            try:
                replaced = replacement(**match_dict)
                if LOG:
                    logger.debug(
                        "%sRule %s.%s: (%s, %s) -> %s", ("  " * (LEVEL)),
                        cls.__name__, key, expr.args, expr.kwargs, replaced)
                return replaced
            except CannotSimplify:
                if LOG_NO_MATCH:
                    logger.debug(
                        "%sRule %s.%s: no match: CannotSimplify",
                        ("  " * (LEVEL)), cls.__name__, key)
                continue
        else:
            if LOG_NO_MATCH:
                logger.debug(
                    "%sRule %s.%s: no match: %s", ("  " * (LEVEL)),
                    cls.__name__, key, match_dict.reason)
    # No matching rules
    return ops, kwargs def match_replace_binary(cls, ops, kwargs):
    """Similar to func:`match_replace`, but for arbitrary length operations,
    such that each two pairs of subsequent operands are matched pairwise.

        >>> A = wc("A")
        >>> class FilterDupes(Operation):
        ...     _binary_rules = {
        ...          'filter_dupes': (pattern_head(A,A), lambda A: A)}
        ...     simplifications = [match_replace_binary, assoc]
        ...     _neutral_element = 0
        >>> FilterDupes.create(1,2,3,4)         # No duplicates
        FilterDupes(1, 2, 3, 4)
        >>> FilterDupes.create(1,2,2,3,4)       # Some duplicates
        FilterDupes(1, 2, 3, 4)

    Note that this only works for *subsequent* duplicate entries:

        >>> FilterDupes.create(1,2,3,2,4)       # No *subsequent* duplicates
        FilterDupes(1, 2, 3, 2, 4)

    Any operation that uses binary reduction must be associative and define a
    neutral element. The binary rules must be compatible with associativity,
    i.e. there is no specific order in which the rules are applied to pairs of
    operands.
    """
    assert assoc in cls.simplifications, (
        cls.__name__ + " must be associative to use match_replace_binary")
    assert hasattr(cls, '_neutral_element'), (
        cls.__name__ + " must define a neutral element to use "
                       "match_replace_binary")
    fops = _match_replace_binary(cls, list(ops))
    if len(fops) == 1:
        return fops[0]
    elif len(fops) == 0:
        return cls._neutral_element
    else:
        return fops, kwargs def orderby(cls, ops, kwargs):
    """Re-order arguments via the class's ``order_key`` key object/function.
    Use this for commutative operations:
    E.g.::

        >>> class Times(Operation):
        ...     order_key = lambda val: val
        ...     simplifications = [orderby, ]
        >>> Times.create(2,1)
        Times(1, 2)
    """
    return sorted(ops, key=cls.order_key), kwargs def delegate_to_method(mtd):
    """Create a simplification rule that delegates the instantiation to the
    method `mtd` of the operand (if defined)"""

    def _delegate_to_method(cls, ops, kwargs):
        assert len(ops) == 1
        op, = ops
        if hasattr(op, mtd):
            return getattr(op, mtd)()
        else:
            return ops, kwargs

    return _delegate_to_method def collect_summands(cls, ops, kwargs):
    """Collect summands that occur multiple times into a single summand

    Also filters out zero-summands.

    Example:
        >>> A, B, C = (OperatorSymbol(s, hs=0) for s in ('A', 'B', 'C'))
        >>> collect_summands(
        ...     OperatorPlus, (A, B, C, ZeroOperator, 2 * A, B, -C) , {})
        ((3 * A^(0), 2 * B^(0)), {})
        >>> collect_summands(OperatorPlus, (A, -A), {})
        ZeroOperator
        >>> collect_summands(OperatorPlus, (B, A, -B), {})
        A^(0)
    """
    from qnet.algebra.core.abstract_quantum_algebra import (
        ScalarTimesQuantumExpression)
    coeff_map = OrderedDict()
    for op in ops:
        if isinstance(op, ScalarTimesQuantumExpression):
            coeff, term = op.coeff, op.term
        else:
            coeff, term = 1, op
        if term in coeff_map:
            coeff_map[term] += coeff
        else:
            coeff_map[term] = coeff
    fops = []
    for (term, coeff) in coeff_map.items():
        op = coeff * term
        if not op.is_zero:
            fops.append(op)
    if len(fops) == 0:
        return cls._zero
    elif len(fops) == 1:
        return fops[0]
    else:
        return tuple(fops), kwargs class BadLiouvillianError(AlgebraError):
    """Raised when a Liouvillian is not of standard Lindblad form."""
    pass class CannotSymbolicallyDiagonalize(AlgebraException):
    """Matrix cannot be diagonalized analytically.

    Signals that a fallback to numerical diagonalization is required.
    """ 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 Matrix(Expression):
    """Matrix of Expressions

    Matrices of :class:`Operator` expressions are required for the SLH
    formalism.
    """
    matrix = None
    _hash = None

    def __init__(self, m):
        if isinstance(m, ndarray):
            self.matrix = m
        elif isinstance(m, Matrix):
            self.matrix = np_array(m.matrix)
        else:
            self.matrix = np_array(m)
        if len(self.matrix.shape) < 2:
            self.matrix = self.matrix.reshape((self.matrix.shape[0], 1))
        if len(self.matrix.shape) > 2:
            raise ValueError()
        super().__init__(self.matrix)

    @property
    def shape(self):
        """The shape of the matrix ``(nrows, ncols)``"""
        return self.matrix.shape

    @property
    def block_structure(self):
        """For square matrices this gives the block (-diagonal) structure of
        the matrix as a tuple of integers that sum up to the full dimension.

        :rtype: tuple
        """
        n, m = self.shape
        if n != m:
            raise AttributeError("block_structure only defined for square "
                                 "matrices")
        for k in range(1, n):
            if ((self.matrix[:k, k:] == 0).all() and
                    (self.matrix[k:, :k] == 0).all()):
                return (k,) + self[k:, k:].block_structure
        return n,

    def _get_blocks(self, block_structure):
        n, m = self.shape
        if n == m:
            if not sum(block_structure) == n:
                raise ValueError()
            if not len(block_structure):
                return ()
            j = block_structure[0]

            if ((self.matrix[:j, j:] == 0).all() and
                    (self.matrix[j:, :j] == 0).all()):
                return ((self[:j, :j],) +
                        self[j:, j:]._get_blocks(block_structure[1:]))
            else:
                raise ValueError()
        elif m == 1:
            if not len(block_structure):
                return ()
            else:
                return ((self[:block_structure[0], :],) +
                        self[:block_structure[0], :]
                        ._get_blocks(block_structure[1:]))
        else:
            raise ValueError()

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

    @property
    def is_zero(self):
        """Are all elements of the matrix zero?"""
        for o in self.matrix.ravel():
            try:
                if not o.is_zero:
                    return False
            except AttributeError:
                if not o == 0:
                    return False
        return True

    @classmethod
    def _get_instance_key(cls, args, kwargs):
        matrix = args[0]
        return (cls, tuple(matrix.ravel()), tuple(matrix.shape))

    def __hash__(self):
        if not self._hash:
            self._hash = hash((tuple(self.matrix.ravel()),
                               self.matrix.shape, Matrix))
        return self._hash

    def __eq__(self, other):
        return (isinstance(other, Matrix) and
                self.shape == other.shape and
                (self.matrix == other.matrix).all())

    def __add__(self, other):
        if isinstance(other, Matrix):
            return Matrix(self.matrix + other.matrix)
        else:
            return Matrix(self.matrix + other)

    def __radd__(self, other):
        return Matrix(other + self.matrix)

    def __mul__(self, other):
        if isinstance(other, Matrix):
            return Matrix(self.matrix.dot(other.matrix))
        else:
            return Matrix(self.matrix * other)

    def __rmul__(self, other):
        return Matrix(other * self.matrix)

    def __sub__(self, other):
        return self + (-1) * other

    def __rsub__(self, other):
        return (-1) * self + other

    def __neg__(self):
        return (-1) * self

    def __truediv__(self, other):
        if is_scalar(other):
            return self * (sympify(1) / other)
        raise NotImplementedError("Can't divide matrix %s by %s"
                                  % (self, other))

    #    def __pow__(self, power):
    #        return OperatorMatrix(self.matrix.__pow__(power))

    def transpose(self):
        """The transpose matrix"""
        return Matrix(self.matrix.T)

    def conjugate(self):
        """The element-wise conjugate matrix

        This is defined only if all the entries in the matrix have a defined
        conjugate (i.e., they have a `conjugate` method). This is *not* the
        case for a matrix of operators. In such a case, only an
        :meth:`elementwise` :func:`adjoint` would be applicable, but this is
        mathematically different from a complex conjugate.

        Raises:
            NoConjugateMatrix: if any entries have no `conjugate` method
        """
        try:
            return Matrix(np_conjugate(self.matrix))
        except AttributeError:
            raise NoConjugateMatrix(
                "Matrix %s contains entries that have no defined "
                "conjugate" % str(self))

    @property
    def real(self):
        """Element-wise real part

        Raises:
            NoConjugateMatrix: if entries have no `conjugate` method and no
                other way to determine the real part

        Note:
            A mathematically equivalent way to obtain a real matrix from a
            complex matrix ``M`` is::

                (M.conjugate() + M) / 2

            However, the result may not be identical to ``M.real``, as the
            latter tries to convert elements of the matrix to real values
            directly, if possible, and only uses the conjugate as a fall-back
        """

        def re(val):
            if hasattr(val, 'real'):
                return val.real
            elif hasattr(val, 'as_real_imag'):
                return val.as_real_imag()[0]
            elif hasattr(val, 'conjugate'):
                return (val.conjugate() + val) / 2
            else:
                raise NoConjugateMatrix(
                    "Matrix entry %s contains has no defined "
                    "conjugate" % str(val))

        # Note: Do NOT use self.matrix.real! This will give wrong results, as
        # numpy thinks of objects (Operators) as real, even if they have no
        # defined real part
        return self.element_wise(re)

    @property
    def imag(self):
        """Element-wise imaginary part

        Raises:
            NoConjugateMatrix: if entries have no `conjugate` method and no
                other way to determine the imaginary part

        Note:
            A mathematically equivalent way to obtain an imaginary matrix from
            a complex matrix ``M`` is::

                (M.conjugate() - M) / (I * 2)

            with same same caveats as :attr:`real`.
        """

        def im(val):
            if hasattr(val, 'imag'):
                return val.imag
            elif hasattr(val, 'as_real_imag'):
                return val.as_real_imag()[1]
            elif hasattr(val, 'conjugate'):
                return (val.conjugate() - val) / (2 * I)
            else:
                raise NoConjugateMatrix(
                    "Matrix entry %s contains has no defined "
                    "conjugate" % str(val))

        # Note: Do NOT use self.matrix.real! This will give wrong results, as
        # numpy thinks of objects (Operators) as real, even if they have no
        # defined real part
        return self.element_wise(im)

    @property
    def T(self):
        """Alias for :meth:`transpose`"""
        return self.transpose()

    def adjoint(self):
        """Adjoint of the matrix

        This is the transpose and the Hermitian adjoint of all elements."""
        return self.T.element_wise(adjoint)

    dag = adjoint

    def trace(self):
        if self.shape[0] == self.shape[1]:
            return sum(self.matrix[k, k] for k in range(self.shape[0]))
        raise NonSquareMatrix(repr(self))

    @property
    def H(self):
        """Alias for :meth:`adjoint`"""
        return self.adjoint()

    def __getitem__(self, item_id):
        item = self.matrix.__getitem__(item_id)
        if isinstance(item, ndarray):
            return Matrix(item)
        return item

    def element_wise(self, func, *args, **kwargs):
        """Apply a function to each matrix element and return the result in a
        new operator matrix of the same shape.

        Args:
            func (FunctionType): A function to be applied to each element. It
                must take the element as its first argument.
            args: Additional positional arguments to be passed to `func`
            kwargs: Additional keyword arguments to be passed to `func`

        Returns:
            Matrix: Matrix with results of `func`, applied element-wise.
        """
        s = self.shape
        emat = [func(o, *args, **kwargs) for o in self.matrix.ravel()]
        return Matrix(np_array(emat).reshape(s))

    def series_expand(self, param: Symbol, about, order: int):
        """Expand the matrix expression as a truncated power series in a scalar
        parameter.

        Args:
            param: Expansion parameter.
            about (.Scalar): Point about which to expand.
            order: Maximum order of expansion >= 0

        Returns:
            tuple of length (order+1), where the entries are the expansion
            coefficients.
        """
        s = self.shape
        emats = zip(*[o.series_expand(param, about, order)
                      for o in self.matrix.ravel()])
        return tuple((Matrix(np_array(em).reshape(s)) for em in emats))

    def expand(self):
        """Expand each matrix element distributively.

        Returns:
            Matrix: Expanded matrix.
        """
        return self.element_wise(
            lambda o: o.expand() if isinstance(o, QuantumExpression) else o)

    def _substitute(self, var_map):
        if self in var_map:
            return var_map[self]
        else:
            return self.element_wise(lambda o: substitute(o, var_map))

    @property
    def free_symbols(self):
        ret = set()
        for o in self.matrix.ravel():
            try:
                ret = ret | o.free_symbols
            except AttributeError:
                pass
        return ret

    @property
    def space(self):
        """Combined Hilbert space of all matrix elements."""
        arg_spaces = [o.space for o in self.matrix.ravel()
                      if hasattr(o, 'space')]
        if len(arg_spaces) == 0:
            return TrivialSpace
        else:
            return ProductSpace.create(*arg_spaces)

    def simplify_scalar(self, func=sympy.simplify):
        """Simplify all scalar expressions appearing in the Matrix."""

        def element_simplify(v):
            if isinstance(v, sympy.Basic):
                return func(v)
            elif isinstance(v, QuantumExpression):
                return v.simplify_scalar(func=func)
            else:
                return v

        return self.element_wise(element_simplify) DENSE_DIMENSION_LIMIT = 1000
II = IdentityOperator
        A = self.A.expand()
        B = self.B.expand()
    A = wc("A", head=Operator)
    B = wc("B", head=Operator)
class Operator(QuantumExpression, metaclass=ABCMeta):
class LocalOperator(Operator, metaclass=ABCMeta):
class OperatorSymbol(QuantumSymbol, Operator):
class IdentityOperator(Operator, metaclass=Singleton):
class ZeroOperator(Operator, metaclass=Singleton):
class LocalSigma(LocalOperator):
class OperatorPlus(QuantumPlus, Operator):
class OperatorTimes(QuantumTimes, Operator):
class ScalarTimesOperator(Operator, ScalarTimesQuantumExpression):
class OperatorDerivative(QuantumDerivative, Operator):
class Commutator(QuantumOperation, Operator):
class OperatorTrace(SingleQuantumOperation, Operator):
class Adjoint(QuantumAdjoint, Operator):
class OperatorPlusMinusCC(SingleQuantumOperation, Operator):
class PseudoInverse(SingleQuantumOperation, Operator):
class NullSpaceProjector(SingleQuantumOperation, Operator):
class OperatorIndexedSum(QuantumIndexedSum, Operator):
    def pseudo_inverse(self):
    def _pseudo_inverse(self):
    def expand_in_basis(self, basis_states=None, hermitian=False):
    def __init__(self, *args, hs):
    def space(self):
    def args(self):
    def kwargs(self):
    def identifier(self):
    def _diff(self, sym):
    def _simplify_scalar(self, func):
    def _pseudo_inverse(self):
    def space(self):
    def args(self):
    def _diff(self, sym):
    def _adjoint(self):
    def _pseudo_inverse(self):
    def space(self):
    def args(self):
    def _diff(self, sym):
    def _adjoint(self):
    def _pseudo_inverse(self):
    def __init__(self, j, k, *, hs):
    def args(self):
    def index_j(self):
    def index_k(self):
    def raise_jk(self, j_incr=0, k_incr=0):
    def _adjoint(self):
    def _pseudo_inverse(self):
def LocalProjector(j, *, hs):
    def _pseudo_inverse(self):
    def _pseudo_inverse(self):
    def has_minus_prefactor(c):
    def _pseudo_inverse(self):
    def __eq__(self, other):
    def __hash__(self):
    def _adjoint(self):
    def _pseudo_inverse(self):
    def __init__(self, A, B):
    def A(self):
    def B(self):
    def doit(self, classes=None, recursive=True, **kwargs):
    def _doit(self, **kwargs):
    def _expand(self):
    def _series_expand(self, param, about, order):
    def _diff(self, sym):
    def _adjoint(self):
    def _pseudo_inverse(self):
    def __init__(self, op, *, over_space):
    def kwargs(self):
    def operand(self):
    def space(self):
    def _expand(self):
    def _series_expand(self, param, about, order):
    def _diff(self, sym):
    def _adjoint(self):
    def _pseudo_inverse(self):
    def _pseudo_inverse(self):
    def __init__(self, op, *, sign=+1):
    def kwargs(self):
    def minimal_kwargs(self):
    def _expand(self):
    def _diff(self, sym):
    def _adjoint(self):
    def _pseudo_inverse(self):
    def doit(self, classes=None, recursive=True, **kwargs):
    def _doit(self, **kwargs):
    def _expand(self):
    def _pseudo_inverse(self):
    def _adjoint(self):
    def _expand(self):
    def _adjoint(self):
    def _pseudo_inverse(self):
    def _pseudo_inverse(self):
def factor_for_trace(ls: HilbertSpace, op: Operator) -> Operator:
def decompose_space(H, A):
def get_coeffs(expr, expand=False, epsilon=0.):
def _coeff_term(op):
def factor_coeff(cls, ops, kwargs):
def adjoint(obj):
def rewrite_with_operator_pm_cc(expr):
    def _combine_operator_p_cc(A, B):
    def _combine_operator_m_cc(A, B):
    def _scal_combine_operator_pm_cc(c, A, d, B): 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 One(Scalar, metaclass=Singleton):
    """The neutral element with respect to scalar multiplication

    Equivalent to the scalar value one::

        >>> One == 1
        True
    """

    _order_name = 'ScalarValue'  # sort like ScalarValue(1)
    _hash_val = 1

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

    @property
    def val(self):
        return self._hash_val

    @property
    def real(self):
        """Real part"""
        return self

    @property
    def imag(self):
        """Imaginary part"""
        return Zero

    @property
    def _order_key(self):
        return KeyTuple([
            self._order_index, self._order_name or self.__class__.__name__,
            self._order_coeff, KeyTuple([self.val, ]), self._order_kwargs])

    def _diff(self, sym):
        return Zero

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

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

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

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

    def _adjoint(self):
        return self

    def __abs__(self):
        return self

    def __neg__(self):
        return ScalarValue(-1)

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

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

    def __mul__(self, other):
        if isinstance(other, self._val_types):
            return ScalarValue.create(other)
        else:
            return other

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

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

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

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

    def __radd__(self, other):
        if isinstance(other, self._val_types):
            return ScalarValue.create(other + 1)
        else:
            return super().__radd__(other)

    def __rsub__(self, other):
        if isinstance(other, self._val_types):
            return ScalarValue.create(other - 1)
        else:
            return super().__radd__(other)

    __rfloordiv = __rtruediv__ = __rmul__ = __rpow__ = __mul__

    def __rmod__(self, other):
        if isinstance(other, self._val_types):
            return ScalarValue.create(other % 1)
        else:
            return super().__rmod__(other)

    def __complex__(self):
        return 1j

    def __int__(self):
        return 1

    def __float__(self):
        return 1.0

    def _sympy_(self):
        return sympy.sympify(1) class DisjunctCommutativeHSOrder():
    """Auxiliary class that generates the correct pseudo-order relation for
    operator products.  Only operators acting on disjoint Hilbert spaces
    are commuted to reflect the order the local factors have in the total
    Hilbert space. I.e., ``sorted(factors, key=DisjunctCommutativeHSOrder)``
    achieves this ordering.
    """

    def __init__(self, op, space_order=None, op_order=None):
        from qnet.algebra.core.hilbert_space_algebra import TrivialSpace
        self.op = op
        self.space = op.space
        if space_order is None:
            self._space_order = op.space._order_key
        else:
            self._space_order = space_order
        if op_order is None:
            self._op_order = expr_order_key(op)
        else:
            self._op_order = op_order
        self.trivial = False
        if op.space is TrivialSpace:
            self.trivial = True

    def __repr__(self):
        from qnet.printing import srepr
        return "%s(%s, space_order=%r, op_order=%r)" % (
                self.__class__.__name__, srepr(self.op), self._space_order,
                self._op_order)

    def __lt__(self, other):
        if self.trivial and other.trivial:
            return self._op_order < other._op_order
        else:
            if self.space.isdisjoint(other.space):
                return self._space_order < other._space_order
        return None  # no ordering class KeyTuple(tuple):
    """A tuple that allows for ordering, facilitating the default ordering of
    Operations. It differs from a normal tuple in that it falls back to string
    comparison if any elements are not directly comparable"""
    def __lt__(self, other):
        for (a, b) in zip(self, other):
            try:
                if a < b:
                    return True
                elif a > b:
                    return False
            except (TypeError, ValueError):
                if str(a) < str(b):
                    return True
                elif str(a) > str(b):
                    return False
        if len(self) < len(other):
            return True
        elif len(self) > len(other):
            return False
        return None

    def __repr__(self):
        return self.__class__.__name__ + tuple.__repr__(self) class Singleton(ABCMeta):
    """Metaclass for singletons

    Any instantiation of a singleton class yields the exact same object, e.g.::

        >>> class MyClass(metaclass=Singleton):
        ...     pass
        >>> a = MyClass()
        >>> b = MyClass()
        >>> a is b
        True

    You can check that an object is a singleton using::

        >>> isinstance(a, SingletonType)
        True
    """
    _instances = {}

    def __call__(cls, *args, **kwargs):
        if cls not in cls._instances:
            cls._instances[cls] = super().__call__(*args, **kwargs)
        return cls._instances[cls]

    @classmethod
    def __instancecheck__(mcs, instance):
        if instance.__class__ is mcs:
            return True
        else:
            return isinstance(instance.__class__, mcs) def singleton_object(cls):
    """Class decorator that transforms (and replaces) a class definition (which
    must have a Singleton metaclass) with the actual singleton object. Ensures
    that the resulting object can still be "instantiated" (i.e., called),
    returning the same object. Also ensures the object can be pickled, is
    hashable, and has the correct string representation (the name of the
    singleton)

    If the class defines a `_hash_val` class attribute, the hash of the
    singleton will be the hash of that value, and the singleton will compare
    equal to that value. Otherwise, the singleton will have a unique hash and
    compare equal only to itself.
    """
    assert isinstance(cls, Singleton), \
        cls.__name__ + " must use Singleton metaclass"

    def self_instantiate(self):
        return self

    cls.__call__ = self_instantiate
    if hasattr(cls, '_hash_val'):
        cls.__hash__ = lambda self: hash(cls._hash_val)
        cls.__eq__ = lambda self, other: other == cls._hash_val
    else:
        cls.__hash__ = lambda self: hash(cls)
        cls.__eq__ = lambda self, other: other is self
    cls.__repr__ = lambda self: cls.__name__
    cls.__reduce__ = lambda self: cls.__name__
    obj = cls()
    obj.__name__ = cls.__name__
    obj.__qualname__ = cls.__qualname__
    return obj
from abc import ABCMeta
from collections import OrderedDict, defaultdict
from numpy import (array as np_array, sqrt as np_sqrt)
from numpy.linalg import eigh
from sympy import (I, Matrix as SympyMatrix, sqrt)
from .abstract_algebra import Operation
from .abstract_quantum_algebra import (
    ScalarTimesQuantumExpression, QuantumExpression, QuantumSymbol,
    QuantumOperation, QuantumPlus, QuantumTimes, QuantumAdjoint,
    QuantumDerivative)
from .algebraic_properties import (
    assoc, filter_neutral, match_replace, match_replace_binary, orderby,
    delegate_to_method, collect_summands)
from .exceptions import BadLiouvillianError, CannotSymbolicallyDiagonalize
from .hilbert_space_algebra import TrivialSpace
from .matrix_algebra import Matrix
from .operator_algebra import (
    Operator, OperatorPlus, ZeroOperator, sympyOne, Adjoint, PseudoInverse)
from .scalar_algebra import is_scalar, One
from ...utils.ordering import DisjunctCommutativeHSOrder, KeyTuple
from ...utils.singleton import Singleton, singleton_object
    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, ]

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


class SuperAdjoint(QuantumAdjoint, SuperOperator):
    r"""Adjoint of a super-operator

    The mathematical notation for this is typically

    .. math::
        {\rm SuperAdjoint}(\mathcal{L}) =: \mathcal{L}^*

    and for any super operator :math:`\mathcal{L}`, its super-adjoint
    :math:`\mathcal{L}^*` satisfies for any pair of operators :math:`M,N`:

    .. math::
        {\rm Tr}[M (\mathcal{L}N)] = Tr[(\mathcal{L}^*M)  N]

    """

    simplifications = [delegate_to_method('_adjoint')]

    def __init__(self, operand):
        super().__init__(operand)


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)