class Expression(metaclass=ABCMeta):
    """Base class for all QNET Expressions

    Expressions should generally be instantiated using the :meth:`create` class
    method, which takes into account the algebraic properties of the Expression
    and and applies simplifications. It also uses memoization to cache all
    known (sub-)expression. This is possible because expressions are intended
    to be immutable. Any changes to an expression should be made through e.g.
    :meth:`substitute` or :meth:`apply_rule`, which returns a new modified
    expression.

    Every expression has a well-defined list of positional and keyword
    arguments that uniquely determine the expression and that may be accessed
    through the :attr:`args` and :attr:`kwargs` property. That is,

    ::

        expr.__class__(*expr.args, **expr.kwargs)

    will return and object identical to `expr`.

    Class attributes:
        instance_caching (bool):  Flag to indicate whether the :meth:`create`
            class method should cache the instantiation of instances. If True,
            repeated calls to :meth:`create` with the same arguments return
            instantly, instead of re-evaluating all simplifications and rules.
        simplifications (list): List of callable simplifications that
            :meth:`create` will use to process its positional and keyword
            arguments. Each callable must take three parameters (the class, the
            list `args` of positional arguments given to :meth:`create` and a
            dictionary `kwargs` of keyword arguments given to :meth:`create`)
            and return either a tuple of new `args` and `kwargs` (which are
            then handed to the next callable), or an :class:`Expression` (which
            is directly returned as the result of the call to :meth:`create`).
            The built-in available simplification callables are in
            :mod:`~qnet.algebra.core.algebraic_properties`
    """
    # Note: all subclasses of Exression that override `__init__` or `create`
    # *must* call the corresponding superclass method *at the end*. Otherwise,
    # caching will not work correctly

    simplifications = []

    # we cache all instances of Expressions for fast construction
    _instances = {}
    instance_caching = True

    # eventually, we should ensure that the create method is idempotent, i.e.
    # expr.create(*expr.args, **expr.kwargs) == expr(*expr.args, **expr.kwargs)
    _create_idempotent = False

    # At this point, match_replace_binary does not yet guarantee this

    def __init__(self, *args, **kwargs):
        self._hash = None
        self._free_symbols = None
        self._bound_symbols = None
        self._all_symbols = None
        self._instance_key = self._get_instance_key(args, kwargs)

    @classmethod
    def create(cls, *args, **kwargs):
        """Instantiate while applying automatic simplifications

        Instead of directly instantiating `cls`, it is recommended to use
        :meth:`create`, which applies simplifications to the args and keyword
        arguments according to the :attr:`simplifications` class attribute, and
        returns an appropriate object (which may or may not be an instance of
        the original `cls`).

        Two simplifications of particular importance are :func:`.match_replace`
        and :func:`.match_replace_binary` which apply rule-based
        simplifications.

        The :func:`.temporary_rules` context manager may be used to allow
        temporary modification of the automatic simplifications that
        :meth:`create` uses, in particular the rules for
        :func:`.match_replace` and :func:`.match_replace_binary`. Inside the
        managed context, the :attr:`simplifications` class attribute may be
        modified and rules can be managed with :meth:`add_rule` and
        :meth:`del_rules`.
        """
        global LEVEL
        if LOG:
            logger = logging.getLogger('QNET.create')
            logger.debug(
                "%s%s.create(*args, **kwargs); args = %s, kwargs = %s",
                ("  " * LEVEL), cls.__name__, args, kwargs)
            LEVEL += 1
        key = cls._get_instance_key(args, kwargs)
        try:
            if cls.instance_caching:
                instance = cls._instances[key]
                if LOG:
                    LEVEL -= 1
                    logger.debug("%s(cached)-> %s", ("  " * LEVEL), instance)
                return instance
        except KeyError:
            pass
        for i, simplification in enumerate(cls.simplifications):
            if LOG:
                try:
                    simpl_name = simplification.__name__
                except AttributeError:
                    simpl_name = "simpl%d" % i
            simplified = simplification(cls, args, kwargs)
            try:
                args, kwargs = simplified
                if LOG:
                    logger.debug(
                        "%s(%s)-> args = %s, kwargs = %s", ("  " * LEVEL),
                        simpl_name, args, kwargs)
            except (TypeError, ValueError):
                # We assume that if the simplification didn't return a tuple,
                # the result is a fully instantiated object
                if cls.instance_caching:
                    cls._instances[key] = simplified
                if cls._create_idempotent and cls.instance_caching:
                    try:
                        key2 = simplified._instance_key
                        if key2 != key:
                            cls._instances[key2] = simplified  # simplified key
                    except AttributeError:
                        #  simplified might e.g. be a scalar and not have
                        #  _instance_key
                        pass
                if LOG:
                    LEVEL -= 1
                    logger.debug(
                        "%s(%s)-> %s", ("  " * LEVEL), simpl_name, simplified)
                return simplified
        if len(kwargs) > 0:
            cls._has_kwargs = True
        instance = cls(*args, **kwargs)
        if cls.instance_caching:
            cls._instances[key] = instance
        if cls._create_idempotent and cls.instance_caching:
            key2 = cls._get_instance_key(args, kwargs)
            if key2 != key:
                cls._instances[key2] = instance  # instantiated key
        if LOG:
            LEVEL -= 1
            logger.debug("%s -> %s", ("  " * LEVEL), instance)
        return instance

    @classmethod
    def _get_instance_key(cls, args, kwargs):
        """Function that calculates a unique "key" (as a tuple) for the
        given args and kwargs. It is the basis of the hash of an Expression,
        and is used for the internal caching of instances. Every Expression
        stores this key in the `_instance_key` attribute.

        Two expressions for which `expr._instance_key` is the same are
        identical by definition (although `expr1 is expr2` generally only holds
        for explicit Singleton instances)
        """
        return (cls,) + nested_tuple(args) + nested_tuple(kwargs)

    @classmethod
    def _rules_attr(cls):
        """Return the name of the attribute with rules for :meth:`create`"""
        from qnet.algebra.core.algebraic_properties import (
            match_replace, match_replace_binary)
        if match_replace in cls.simplifications:
            return '_rules'
        elif match_replace_binary in cls.simplifications:
            return '_binary_rules'
        else:
            raise TypeError(
                "class %s does not have match_replace or "
                "match_replace_binary in its simplifications" % cls.__name__)

    @classmethod
    def add_rule(cls, name, pattern, replacement, attr=None):
        """Add an algebraic rule for :meth:`create` to the class

        Args:
            name (str): Name of the rule. This is used for debug logging to
                allow an analysis of which rules where applied when creating an
                expression. The `name` can be arbitrary, but it must be unique.
                Built-in rules have names ``'Rxxx'`` where ``x`` is a digit
            pattern (.Pattern): A pattern constructed by :func:`.pattern_head`
                to match a :class:`.ProtoExpr`
            replacement (callable): callable that takes the wildcard names
                defined in `pattern` as keyword arguments and returns an
                evaluated expression.
            attr (None or str): Name of the class attribute to which to add the
                rule. If None, one of ``'_rules'``, ``'_binary_rules'`` is
                automatically chosen

        Raises:
            TypeError: if `name` is not a :class:`str` or `pattern` is not a
                :class:`.Pattern` instance
            ValueError: if `pattern` is not set up to match a
                :class:`.ProtoExpr`; if there there is already a rule with the
                same `name`; if `replacement` is not a callable or does not
                take all the wildcard names in `pattern` as arguments
            AttributeError: If invalid `attr`

        Note:
            The "automatic" rules added by this method are applied *before*
            expressions are instantiated (against a corresponding
            :class:`.ProtoExpr`). In contrast,
            :meth:`apply_rules`/:meth:`apply_rule` are applied to fully
            instantiated objects.

            The :func:`.temporary_rules` context manager may be used to create
            a context in which rules may be defined locally.
        """
        from qnet.utils.check_rules import check_rules_dict
        if attr is None:
            attr = cls._rules_attr()
        if name in getattr(cls, attr):
            raise ValueError(
                "Duplicate key '%s': rule already exists" % name)
        getattr(cls, attr).update(check_rules_dict(
            [(name, (pattern, replacement))]))

    @classmethod
    def show_rules(cls, *names, attr=None):
        """Print algebraic rules used by :class:`create`

        Print a summary of the algebraic rules with the given names, or all
        rules if not names a given.

        Args:
            names (str): Names of rules to show
            attr (None or str): Name of the class attribute from which to get
                the rules. Cf. :meth:`add_rule`.

        Raises:
            AttributeError: If invalid `attr`
        """
        from qnet.printing import srepr
        try:
            if attr is None:
                attr = cls._rules_attr()
            rules = getattr(cls, attr)
        except TypeError:
            rules = {}
        for (name, rule) in rules.items():
            if len(names) > 0 and name not in names:
                continue
            pat, repl = rule
            print(name)
            print("    PATTERN:")
            print(textwrap.indent(
                textwrap.dedent(srepr(pat, indented=True)),
                prefix=" "*8))
            print("    REPLACEMENT:")
            print(textwrap.indent(
                textwrap.dedent(inspect.getsource(repl).rstrip()),
                prefix=" "*8))

    @classmethod
    def del_rules(cls, *names, attr=None):
        """Delete algebraic rules used by :meth:`create`

        Remove the rules with the given `names`, or all rules if no names are
        given

        Args:
            names (str): Names of rules to delete
            attr (None or str): Name of the class attribute from which to
                delete the rules. Cf. :meth:`add_rule`.

        Raises:
            KeyError: If any rules in `names` does not exist
            AttributeError: If invalid `attr`
        """
        if attr is None:
            attr = cls._rules_attr()
        if len(names) == 0:
            getattr(cls, attr)  # raise AttributeError if wrong attr
            setattr(cls, attr, OrderedDict())
        else:
            for name in names:
                del getattr(cls, attr)[name]

    @classmethod
    def rules(cls, attr=None):
        """Iterable of rule names used by :meth:`create`

        Args:
            attr (None or str): Name of the class attribute to which to get the
                names. If None, one of ``'_rules'``, ``'_binary_rules'`` is
                automatically chosen
        """
        try:
            if attr is None:
                attr = cls._rules_attr()
            return getattr(cls, attr).keys()
        except TypeError:
            return ()

    @property
    @abstractmethod
    def args(self):
        """The tuple of positional arguments for the instantiation of the
        Expression"""
        raise NotImplementedError(self.__class__.__name__)

    @property
    def kwargs(self):
        """The dictionary of keyword-only arguments for the instantiation of
        the Expression"""
        # Subclasses must override this property if and only if they define
        # keyword-only arguments in their __init__ method
        if hasattr(self, '_has_kwargs') and self._has_kwargs:
            raise NotImplementedError(
                "Class %s does not provide a kwargs property"
                % str(self.__class__.__name__))
        return {}

    @property
    def minimal_kwargs(self):
        """A "minimal" dictionary of keyword-only arguments, i.e. a subset of
        `kwargs` that may exclude default options"""
        return self.kwargs

    def __eq__(self, other):
        try:
            return ((self is other) or
                    (self._instance_key == other._instance_key))
        except AttributeError:
            return False

    def __hash__(self):
        if self._hash is None:
            self._hash = hash(self._instance_key)
        return self._hash

    def __repr__(self):
        # This method will be replaced by init_printing()
        from qnet.printing import init_printing
        init_printing()
        return repr(self)

    def __str__(self):
        # This method will be replaced by init_printing()
        from qnet.printing import init_printing
        init_printing()
        return str(self)

    def substitute(self, var_map):
        """Substitute sub-expressions

        Args:
            var_map (dict): Dictionary with entries of the form
                ``{expr: substitution}``
        """
        if self in var_map:
            return var_map[self]
        return self._substitute(var_map)

    def _substitute(self, var_map, safe=False):
        """Implementation of :meth:`substitute`.

        For internal use, the `safe` keyword argument allows to perform a
        substitution on the `args` and `kwargs` of the expression only,
        guaranteeing that the type of the expression does not change, at the
        cost of possibly not returning a maximally simplified expression. The
        `safe` keyword is not handled recursively, i.e. any `args`/`kwargs`
        will be fully simplified, possibly changing their types.
        """
        if self in var_map:
            if not safe or (type(var_map[self]) == type(self)):
                return var_map[self]
        if isinstance(self.__class__, Singleton):
            return self
        new_args = [substitute(arg, var_map) for arg in self.args]
        new_kwargs = {key: substitute(val, var_map)
                      for (key, val) in self.kwargs.items()}
        if safe:
            return self.__class__(*new_args, **new_kwargs)
        else:
            return self.create(*new_args, **new_kwargs)

    def doit(self, classes=None, recursive=True, **kwargs):
        """Rewrite (sub-)expressions in a more explicit form

        Return a modified expression that is more explicit than the original
        expression. The definition of "more explicit" is decided by the
        relevant subclass, e.g. a :meth:`Commutator <.Commutator.doit>` is
        written out according to its definition.

        Args:
            classes (None or list): an optional list of classes. If given,
                only (sub-)expressions that an instance of one of the classes
                in the list will be rewritten.
            recursive (bool): If True, also rewrite any sub-expressions of any
                rewritten expression. Note that :meth:`doit` always recurses
                into sub-expressions of expressions not affected by it.
            kwargs: Any remaining keyword arguments may be used by the
                :meth:`doit` method of a particular expression.

        Example:

            Consider the following expression::

                >>> from sympy import IndexedBase
                >>> i = IdxSym('i'); N = symbols('N')
                >>> Asym, Csym = symbols('A, C', cls=IndexedBase)
                >>> A = lambda i: OperatorSymbol(StrLabel(Asym[i]), hs=0)
                >>> B = OperatorSymbol('B', hs=0)
                >>> C = lambda i: OperatorSymbol(StrLabel(Csym[i]), hs=0)
                >>> def show(expr):
                ...     print(unicode(expr, show_hs_label=False))
                >>> expr = Sum(i, 1, 3)(Commutator(A(i), B) + C(i)) / N
                >>> show(expr)
                1/N (∑_{i=1}^{3} (Ĉ_i + [Â_i, B̂]))

            Calling :meth:`doit` without parameters rewrites both the indexed
            sum and the commutator::

                >>> show(expr.doit())
                1/N (Ĉ₁ + Ĉ₂ + Ĉ₃ + Â₁ B̂ + Â₂ B̂ + Â₃ B̂ - B̂ Â₁ - B̂ Â₂ - B̂ Â₃)

            A non-recursive call only expands the sum, as it does not recurse
            into the expanded summands::

                >>> show(expr.doit(recursive=False))
                1/N (Ĉ₁ + Ĉ₂ + Ĉ₃ + [Â₁, B̂] + [Â₂, B̂] + [Â₃, B̂])

            We can selectively expand only the sum or only the commutator::

                >>> show(expr.doit(classes=[IndexedSum]))
                1/N (Ĉ₁ + Ĉ₂ + Ĉ₃ + [Â₁, B̂] + [Â₂, B̂] + [Â₃, B̂])

                >>> show(expr.doit(classes=[Commutator]))
                1/N (∑_{i=1}^{3} (Ĉ_i - B̂ Â_i + Â_i B̂))

            Also we can pass a keyword argument that expands the sum only to
            the 2nd term, as documented in :meth:`.Commutator.doit`

                >>> show(expr.doit(classes=[IndexedSum], max_terms=2))
                1/N (Ĉ₁ + Ĉ₂ + [Â₁, B̂] + [Â₂, B̂])
        """
        in_classes = (
            (classes is None) or
            any([isinstance(self, cls) for cls in classes]))
        if in_classes:
            new = self._doit(**kwargs)
        else:
            new = self
        if (new == self) or recursive:
            new_args = []
            for arg in new.args:
                if isinstance(arg, Expression):
                    new_args.append(arg.doit(
                        classes=classes, recursive=recursive, **kwargs))
                else:
                    new_args.append(arg)
            new_kwargs = OrderedDict([])
            for (key, val) in new.kwargs.items():
                if isinstance(val, Expression):
                    new_kwargs[key] = val.doit(
                        classes=classes, recursive=recursive, **kwargs)
                else:
                    new_kwargs[key] = val
            new = new.__class__.create(*new_args, **new_kwargs)
            if new != self and recursive:
                new = new.doit(classes=classes, recursive=True, **kwargs)
        return new

    def _doit(self, **kwargs):
        """Non-recursively rewrite expression in a more explicit form"""
        # Any subclass that overrides :meth:`_doit` should also override
        # :meth:`doit` with a stub (calling ``super().doit`` only), but
        # also provide the documentation for :meth:`_doit` (since :meth:`_doit`
        # won't be rendered by Sphinx)
        return self

    def apply(self, func, *args, **kwargs):
        """Apply `func` to expression.

        Equivalent to ``func(self, *args, **kwargs)``. This method exists for
        easy chaining::

            >>> A, B, C, D = (
            ...     OperatorSymbol(s, hs=1) for s in ('A', 'B', 'C', 'D'))
            >>> expr = (
            ...     Commutator(A * B, C * D)
            ...     .apply(lambda expr: expr**2)
            ...     .apply(expand_commutators_leibniz, expand_expr=False)
            ...     .substitute({A: IdentityOperator}))
        """
        return func(self, *args, **kwargs)

    def apply_rules(self, rules, recursive=True):
        """Rebuild the expression while applying a list of rules

        The rules are applied against the instantiated expression, and any
        sub-expressions if `recursive` is True. Rule application is best though
        of as a pattern-based substitution. This is different from the
        *automatic* rules that :meth:`create` uses (see :meth:`add_rule`),
        which are applied *before* expressions are instantiated.

        Args:
            rules (list or ~collections.OrderedDict): List of rules or
                dictionary mapping names to rules, where each rule is a tuple
                (:class:`Pattern`, replacement callable), cf.
                :meth:`apply_rule`
            recursive (bool): If true (default), apply rules to all arguments
                and keyword arguments of the expression. Otherwise, only the
                expression itself will be re-instantiated.

        If `rules` is a dictionary, the keys (rules names) are used only for
        debug logging, to allow an analysis of which rules lead to the final
        form of an expression.
        """
        if recursive:
            new_args = [_apply_rules(arg, rules) for arg in self.args]
            new_kwargs = {
                key: _apply_rules(val, rules)
                for (key, val) in self.kwargs.items()}
        else:
            new_args = self.args
            new_kwargs = self.kwargs
        simplified = self.create(*new_args, **new_kwargs)
        return _apply_rules_no_recurse(simplified, rules)

    def apply_rule(self, pattern, replacement, recursive=True):
        """Apply a single rules to the expression

        This is equivalent to :meth:`apply_rules` with
        ``rules=[(pattern, replacement)]``

        Args:
            pattern (.Pattern): A pattern containing one or more wildcards
            replacement (callable): A callable that takes the wildcard names in
                `pattern` as keyword arguments, and returns a replacement for
                any expression that `pattern` matches.

        Example:
            Consider the following Heisenberg Hamiltonian::

                >>> tls = SpinSpace(label='s', spin='1/2')
                >>> i, j, n = symbols('i, j, n', cls=IdxSym)
                >>> J = symbols('J', cls=sympy.IndexedBase)
                >>> def Sig(i):
                ...     return OperatorSymbol(
                ...         StrLabel(sympy.Indexed('sigma', i)), hs=tls)
                >>> H = - Sum(i, tls)(Sum(j, tls)(
                ...     J[i, j] * Sig(i) * Sig(j)))
                >>> unicode(H)
                '- (∑_{i,j ∈ ℌₛ} J_ij σ̂_i^(s) σ̂_j^(s))'

            We can transform this into a classical Hamiltonian by replacing the
            operators with scalars::

                >>> H_classical = H.apply_rule(
                ...     pattern(OperatorSymbol, wc('label', head=StrLabel)),
                ...     lambda label: label.expr * IdentityOperator)
                >>> unicode(H_classical)
                '- (∑_{i,j ∈ ℌₛ} J_ij σ_i σ_j)'
        """
        return self.apply_rules([(pattern, replacement)], recursive=recursive)

    def rebuild(self):
        """Recursively re-instantiate the expression

        This is generally used within a managed context such as
        :func:`.extra_rules`, :func:`.extra_binary_rules`, or
        :func:`.no_rules`.
        """
        return self.apply_rules(rules={})

    def _repr_latex_(self):
        """For compatibility with the IPython notebook, generate TeX expression
        and surround it with $'s.
        """
        # This method will be replaced by init_printing()
        from qnet.printing import init_printing
        init_printing()
        return self._repr_latex_()

    def _sympy_(self):
        # By default, when a QNET expression occurring in a SymPy context (e.g.
        # when converting a QNET Matrix to a Sympy Matrix), sympify will try to
        # parse the string representation of the Expression. This will usually
        # fail, but when it doesn't, it always produces nonsense. Thus, we make
        # it fail explicitly
        raise SympifyError("QNET expressions cannot be converted to SymPy")

    @property
    def free_symbols(self):
        """Set of free SymPy symbols contained within the expression."""
        if self._free_symbols is None:
            res = set.union(
                set([]), # dummy arg (union fails without arguments)
                *[_free_symbols(val) for val in self.kwargs.values()])
            res.update(
                set([]), # dummy arg (update fails without arguments)
                *[_free_symbols(arg) for arg in self.args])
            self._free_symbols = res
        return self._free_symbols

    @property
    def bound_symbols(self):
        """Set of bound SymPy symbols in the expression"""
        if self._bound_symbols is None:
            res = set.union(
                set([]), # dummy arg (union fails without arguments)
                *[_bound_symbols(val) for val in self.kwargs.values()])
            res.update(
                set([]), # dummy arg (update fails without arguments)
                *[_bound_symbols(arg) for arg in self.args])
            self._bound_symbols = res
        return self._bound_symbols

    @property
    def all_symbols(self):
        """Combination of :attr:`free_symbols` and :attr:`bound_symbols`"""
        if self._all_symbols is None:
            self._all_symbols = self.free_symbols | self.bound_symbols
        return self._all_symbols

    def __ne__(self, other):
        """If it is well-defined (i.e. boolean), simply return
        the negation of ``self.__eq__(other)``
        Otherwise return NotImplemented.
        """
        eq = self.__eq__(other)
        if type(eq) is bool:
            return not eq
        return NotImplemented 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 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 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 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 QuantumIndexedSum(IndexedSum, SingleQuantumOperation, metaclass=ABCMeta):
    """Base class for indexed sums"""

    @property
    def space(self):
        """The Hilbert space of the sum's term"""
        return self.term.space

    def _expand(self):
        return self.__class__.create(self.term.expand(), *self.ranges)

    def _adjoint(self):
        return self.__class__.create(self.term.adjoint(), *self.ranges)

    def _diff(self, sym):
        return self.__class__.create(self.term.diff(sym), *self.ranges)

    def __mul__(self, other):
        from qnet.algebra.core.scalar_algebra import is_scalar
        if isinstance(other, IndexedSum):
            other = other.make_disjunct_indices(self)
            new_ranges = self.ranges + other.ranges
            new_term = self.term * other.term
            # note that class may change, depending on type of new_term
            return new_term.__class__._indexed_sum_cls.create(
                new_term, *new_ranges)
        elif is_scalar(other):
            return self.__class__._scalar_times_expr_cls(other, self)
        elif isinstance(other, ScalarTimesQuantumExpression):
            return self._class__._scalar_times_expr_cls(
                other.coeff, self * other.term)
        else:
            sum = self.make_disjunct_indices(*other.bound_symbols)
            new_term = sum.term * other
            return new_term.__class__._indexed_sum_cls.create(
                new_term, *sum.ranges)

    def __rmul__(self, other):
        from qnet.algebra.core.scalar_algebra import is_scalar
        if isinstance(other, IndexedSum):
            self_new = self.make_disjunct_indices(other)
            new_ranges = other.ranges + self_new.ranges
            new_term = other.term * self_new.term
            # note that class may change, depending on type of new_term
            return new_term.__class__._indexed_sum_cls.create(
                new_term, *new_ranges)
        elif is_scalar(other):
            return self.__class__._scalar_times_expr_cls.create(other, self)
        elif isinstance(other, ScalarTimesQuantumExpression):
            return self._class__._scalar_times_expr_cls(
                other.coeff, other.term * self)
        else:
            sum = self.make_disjunct_indices(*other.bound_symbols)
            new_term = other * sum.term
            return new_term.__class__._indexed_sum_cls.create(
                new_term, *sum.ranges)

    def __add__(self, other):
        raise NotImplementedError()

    def __radd__(self, other):
        raise NotImplementedError()

    def __sub__(self, other):
        raise NotImplementedError()

    def __rsub__(self, other):
        raise NotImplementedError() 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 ensure_local_space(hs, cls=LocalSpace):
    """Ensure that the given `hs` is an instance of :class:`LocalSpace`.

    If `hs` an instance of :class:`str` or :class:`int`, it will be converted
    to a `cls` (if possible). If it already is an instace of `cls`, `hs`
    will be returned unchanged.

    Args:
        hs (HilbertSpace or str or int): The Hilbert space (or label) to
            convert/check
        cls (type): The class to which an int/str label for a Hilbert space
            should be converted. Must be a subclass of :class:`LocalSpace`.

    Raises:
        TypeError: If `hs` is not a :class:`.LocalSpace`, :class:`str`, or
            :class:`int`.

    Returns:
        LocalSpace: original or converted `hs`

    Examples:
        >>> srepr(ensure_local_space(0))
        "LocalSpace('0')"
        >>> srepr(ensure_local_space('tls'))
        "LocalSpace('tls')"
        >>> srepr(ensure_local_space(0, cls=LocalSpace))
        "LocalSpace('0')"
        >>> srepr(ensure_local_space(LocalSpace(0)))
        "LocalSpace('0')"
        >>> srepr(ensure_local_space(LocalSpace(0)))
        "LocalSpace('0')"
        >>> srepr(ensure_local_space(LocalSpace(0) * LocalSpace(1)))
        Traceback (most recent call last):
           ...
        TypeError: hs must be an instance of LocalSpace
    """
    if isinstance(hs, (str, int)):
        try:
            hs = cls(hs)
        except TypeError as exc_info:
            raise TypeError(
                "Cannot convert %s '%s' into a %s instance: %s"
                % (hs.__class__.__name__, hs, cls.__name__, str(exc_info)))
    if not isinstance(hs, LocalSpace):
        raise TypeError("hs must be an instance of LocalSpace")
    return hs def _series_expand_combine_prod(c1, c2, order):
    """Given the result of the ``c1._series_expand(...)`` and
    ``c2._series_expand(...)``, construct the result of
    ``(c1*c2)._series_expand(...)``
    """
    from qnet.algebra.core.scalar_algebra import Zero
    res = []
    c1 = list(c1)
    c2 = list(c2)
    for n in range(order + 1):
        summands = []
        for k in range(n + 1):
            if c1[k].is_zero or c2[n-k].is_zero:
                summands.append(Zero)
            else:
                summands.append(c1[k] * c2[n - k])
        sum = summands[0]
        for summand in summands[1:]:
            if summand != 0:
                sum += summand
        res.append(sum)
    return tuple(res) def accept_bras(cls, ops, kwargs):
    """Accept operands that are all bras, and turn that into to bra of the
    operation applied to all corresponding kets"""
    from qnet.algebra.core.state_algebra import Bra
    kets = []
    for bra in ops:
        if isinstance(bra, Bra):
            kets.append(bra.ket)
        else:
            return ops, kwargs
    return Bra.create(cls.create(*kets, **kwargs)) 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 assoc_indexed(cls, ops, kwargs):
    r"""Flatten nested indexed structures while pulling out possible prefactors

    For example, for an :class:`.IndexedSum`:

    .. math::

        \sum_j \left( a \sum_i \dots \right) = a \sum_{j, i} \dots
    """
    from qnet.algebra.core.abstract_quantum_algebra import (
        ScalarTimesQuantumExpression)
    term, *ranges = ops

    if isinstance(term, cls):
        coeff = 1
    elif isinstance(term, ScalarTimesQuantumExpression):
        coeff = term.coeff
        term = term.term
        if not isinstance(term, cls):
            return ops, kwargs
    else:
        return ops, kwargs

    term = term.make_disjunct_indices(*ranges)
    combined_ranges = tuple(ranges) + term.ranges

    if coeff == 1:
        return cls.create(term.term, *combined_ranges)
    else:
        bound_symbols = set([r.index_symbol for r in combined_ranges])
        if len(coeff.free_symbols.intersection(bound_symbols)) == 0:
            return coeff * cls.create(term.term, *combined_ranges)
        else:
            return cls.create(coeff * term.term, *combined_ranges) def basis_ket_zero_outside_hs(cls, ops, kwargs):
    """For ``BasisKet.create(ind, hs)`` with an integer label `ind`, return a
    :class:`ZeroKet` if `ind` is outside of the range of the underlying Hilbert
    space
    """
    from qnet.algebra.core.state_algebra import ZeroKet
    ind, = ops
    hs = kwargs['hs']
    if isinstance(ind, int):
        if ind < 0 or (hs._dimension is not None and ind >= hs._dimension):
            return ZeroKet
    return ops, 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 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 OverlappingSpaces(AlgebraError):
    """Raised when objects fail to be in separate Hilbert spaces.""" class UnequalSpaces(AlgebraError):
    """Raised when objects fail to be in the same Hilbert space.

    This happens for example when trying to add two states from different
    Hilbert spaces.""" class SpaceTooLargeError(AlgebraError):
    """Raised when objects fail to be have overlapping Hilbert spaces.""" class FullSpace(HilbertSpace, Expression, metaclass=Singleton):
    """The 'full space', i.e. a Hilbert space that includes any other Hilbert
    space as a tensor factor.

    The `FullSpace` has no defined basis, any related properties will raise
    :class:`.BasisNotSetError`
    """

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

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

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

    def remove(self, other):
        """Raise AlgebraError, as the remaining space is undefined"""
        raise AlgebraError("Cannot remove anything from FullSpace")

    @property
    def local_factors(self):
        """Raise AlgebraError, as the the factors of the full space are
        undefined"""
        raise AlgebraError("FullSpace has no local_factors")

    def intersect(self, other):
        """Return `other`"""
        return other

    def is_strict_subfactor_of(self, other):
        """False, as the full space by definition is not contained in any other
        space"""
        return False

    @property
    def label(self):
        return "total" 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 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 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 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 ScalarExpression(Scalar, metaclass=ABCMeta):
    """Base class for scalars with non-scalar arguments

    For example, a :class:`.BraKet` is a :class:`Scalar`, but has arguments
    that are states.
    """

    _order_index = 1  # Expression scalars come after ScalarValue

    def __pow__(self, other):
        return ScalarPower.create(self, other) class FockIndex(IntIndex):
    """Symbolic index labeling a basis state in a :class:`.LocalSpace`"""

    @property
    def fock_index(self):
        return self.expr class IdxSym(sympy.Symbol):
    """Index symbol in an indexed sum or product

    Args:
        name (str): The label for the symbol. It must be a simple Latin or
            Greek letter, possibly with a subscript, e.g. ``'i'``, ``'mu'``,
            ``'gamma_A'``
        primed (int): Number of prime marks (') associated with the symbol

    Notes:

        The symbol can be used in arbitrary algebraic (sympy) expressions::

            >>> sympy.sqrt(IdxSym('n') + 1)
            sqrt(n + 1)

        By default, the symbol is assumed to represent an integer. If this is
        not the case, you can instantiate explicitly as a non-integer::

            >>> IdxSym('i').is_integer
            True
            >>> IdxSym('i', integer=False).is_integer
            False

        You may also declare the symbol as positive::

            >>> IdxSym('i').is_positive
            >>> IdxSym('i', positive=True).is_positive
            True

        The `primed` parameter is used to automatically create distinguishable
        indices in products of sums, or more generally if the same index occurs
        in an expression with potentially differnt values::

            >>> ascii(IdxSym('i', primed=2))
            "i''"
            >>> IdxSym('i') == IdxSym('i', primed=1)
            False

        It should not be used when creating indices "by hand"

    Raises:
        ValueError: if `name` is not a simple symbol label, or if primed < 0
        TypeError: if `name` is not a string
    """

    is_finite = True
    is_Symbol = True
    is_symbol = True
    is_Atom = True
    _diff_wrt = True

    _rx_name = re.compile('^[A-Za-z]+(_[A-Za-z0-9().,+-]+)?$')

    def __new_stage2__(cls, name, primed=0, **kwargs):
        # remove: start, stop, points
        if not cls._rx_name.match(name):
            raise ValueError(
                "name '%s' does not match pattern '%s'"
                % (name, cls._rx_name.pattern))
        primed = int(primed)
        if not primed >= 0:
            raise ValueError("`primed` must be an integer >= 0")
        if 'integer' not in kwargs:
            kwargs['integer'] = True
        obj = super().__xnew__(cls, name, **kwargs)
        obj.params = (primed, )
        obj._primed = primed
        return obj

    def __new__(cls, name, *, primed=0, **kwargs):
        obj = cls.__xnew_cached_(cls, name, primed, **kwargs)
        return obj

    __xnew__ = staticmethod(__new_stage2__)
    __xnew_cached_ = staticmethod(sympy_cacheit(__new_stage2__))

    def _hashable_content(self):
        return (sympy.Symbol._hashable_content(self), self.params)

    @property
    def primed(self):
        return self._primed

    def incr_primed(self, incr=1):
        """Return a copy of the index with an incremented :attr:`primed`"""
        return self.__class__(
            self.name, primed=self._primed + incr,
            **self._assumptions.generator)

    @property
    def prime(self):
        """equivalent to :meth:`inc_primed` with ``incr=1``"""
        return self.incr_primed(incr=1)

    def _sympystr(self, printer, *args):
        return printer._print_Symbol(self) + "'" * self.primed

    def _sympyrepr(self, printer, *args):
        res = printer._print_Symbol(self)
        if self.primed > 0:
            res = res[:-1] + ", primed=%d)" % self.primed
        return res

    def _pretty(self, printer, *args):
        res = printer._print_Symbol(self)
        return res.__class__(
            res.s + "'" * self.primed, res.baseline, res.binding)

    def _latex(self, printer, *args):
        res = printer._print_Symbol(self)
        if self.primed > 0:
            res = r'{%s^{%s}}' % (res, r'\prime' * self.primed)
        return res class IndexOverFockSpace(IndexRangeBase):
    """Index range over the integer indices of a :class:`.LocalSpace` basis

    Args:
        index_symbol (IdxSym): The symbol iterating over the range
        hs (.LocalSpace): Hilbert space over whose basis to iterate

    """
    hs = attr.ib()
    # TODO: assert that hs is indeed a LocalSpace

    def iter(self):
        if self.hs._dimension is None:
            i = 0
            while True:
                yield {self.index_symbol: i}
                i += 1
        else:
            for ind in range(self.hs.dimension):
                yield {self.index_symbol: ind}

    def __len__(self):
        return self.hs._dimension

    def __contains__(self, val):
        if self.hs._dimension is None:
            return val >= 0
        else:
            return 0 <= val < self.hs.dimension

    def substitute(self, var_map):
        new_index_symbol = var_map.get(self.index_symbol, self.index_symbol)
        new_hs = var_map.get(self.hs, self.hs)
        return self.__class__(index_symbol=new_index_symbol, hs=new_hs)

    def piecewise_one(self, expr):
        if self.hs._dimension is None:
            to_val = sympy.oo
        else:
            to_val = self.hs.dimension
        return Piecewise(
            (1, sympy.Interval(0, to_val).as_relational(expr)),
            (0, True)) class IndexOverRange(IndexRangeBase):
    """Index over the inclusive range between two integers

    Args:
        index_symbol (IdxSym): The symbol iterating over the range
        start_from (int): Starting value for the index
        to (int): End value of the index
        step (int): Step width by which index increases
    """
    start_from = attr.ib(validator=attr.validators.instance_of(int))
    to = attr.ib(validator=attr.validators.instance_of(int))
    step = attr.ib(validator=attr.validators.instance_of(int), default=1)

    def iter(self):
        for ind in self.range:
            yield {self.index_symbol: ind}

    @property
    def range(self):
        return range(
            self.start_from,
            (self.to + 1) if self.step >= 0 else (self.to - 1),
            self.step)

    def __len__(self):
        return len(self.range)

    def __contains__(self, val):
        return val in self.range

    def substitute(self, var_map):
        new_index_symbol = var_map.get(self.index_symbol, self.index_symbol)
        return self.__class__(
            index_symbol=new_index_symbol, start_from=self.start_from,
            to=self.to, step=self.step)

    def piecewise_one(self, expr):
        return Piecewise(
            (1, sympy.FiniteSet(*list(self.range)).contains(expr)),
            (0, True)) class SymbolicLabelBase(metaclass=ABCMeta):
    """Base class for symbolic labels

    A symbolic label is a SymPy expression that contains one or more
    :class:`IdxSym`, and can be rendered into an integer or string label by
    substituting integer values for each :class:`IdxSym`.

    See :class:`IntIndex` for an example.
    """
    expr = attr.ib()

    @expr.validator
    def _validate_expr(self, attribute, value):
        if not isinstance(value, sympy.Basic):
            raise ValueError("expr must be a sympy formula")
        if not self._has_idx_syms(value):
            raise ValueError("expr must contain at least one IdxSym")

    @staticmethod
    def _has_idx_syms(expr):
        return any([isinstance(sym, IdxSym) for sym in expr.free_symbols])

    @abstractmethod
    def _render(self, expr):
        """Render `expr` into a a label. It can be assumed that `expr` does not
        contain any :class:`IdxSym`"""
        pass

    def substitute(self, var_map):
        """Substitute in the expression describing the label.

        If the result of the substitution no longer contains any
        :class:`IdxSym`, this returns a "rendered" label.
        """
        new_expr = self.expr.subs(var_map)
        if self._has_idx_syms(new_expr):
            return self.__class__(expr=new_expr)
        else:
            return self._render(new_expr)

    @property
    def free_symbols(self):
        """Free symbols in the expression describing the label"""
        return self.expr.free_symbols class FullCommutativeHSOrder():
    """Auxiliary class that generates the correct pseudo-order relation for
    operator sums.  Operators are first ordered by their Hilbert space, then by
    their order-key; ``sorted(factors, key=FullCommutativeHSOrder)``
    achieves this ordering.
    """

    def __init__(self, op, space_order=None, op_order=None):
        self.op = op
        self.space = op.space
        if space_order is None:
            self._space_order = self.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

    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.space == other.space:
            return self._op_order < other._op_order
        else:
            return self._space_order < other._space_order 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
import re
import sympy
from abc import ABCMeta
from collections import OrderedDict
from itertools import product as cartesian_product
from .abstract_algebra import Expression, Operation
from .scalar_algebra import is_scalar
from .abstract_quantum_algebra import (
    ScalarTimesQuantumExpression, QuantumExpression, QuantumSymbol,
    QuantumPlus, QuantumTimes, QuantumAdjoint, QuantumIndexedSum,
    QuantumDerivative, ensure_local_space, _series_expand_combine_prod)
from .algebraic_properties import (
    accept_bras, assoc, assoc_indexed, basis_ket_zero_outside_hs,
    filter_neutral, match_replace, match_replace_binary, orderby,
    collect_summands)
from .exceptions import OverlappingSpaces, UnequalSpaces, SpaceTooLargeError
from .hilbert_space_algebra import FullSpace, TrivialSpace
from qnet.algebra.core.algebraic_properties import (
    indexed_sum_over_const,
    indexed_sum_over_kronecker)
from .operator_algebra import Operator, OperatorPlus, PseudoInverse
from .scalar_algebra import ScalarExpression
from ...utils.indices import (
    FockIndex, IdxSym, IndexOverFockSpace, IndexOverRange, SymbolicLabelBase, )
from ...utils.ordering import FullCommutativeHSOrder
from ...utils.singleton import Singleton, singleton_object
        from qnet.algebra.library.fock_operators import Destroy, Create
    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):
        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)


###############################################################################
# Algebra Operations
###############################################################################


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]