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 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 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 check_cdims(cls, ops, kwargs):
    """Check that all operands (`ops`) have equal channel dimension."""
    if not len({o.cdim for o in ops}) == 1:
        raise ValueError("Not all operands have the same cdim:" + str(ops))
    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 filter_cid(cls, ops, kwargs):
    """Remove occurrences of the :func:`.circuit_identity` ``cid(n)`` for any
    ``n``. Cf. :func:`filter_neutral`
    """
    from qnet.algebra.core.circuit_algebra import CircuitZero, circuit_identity
    if len(ops) == 0:
        return CircuitZero
    fops = [op for op in ops if op != circuit_identity(op.cdim)]
    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 class AlgebraError(AlgebraException):
    """Base class for all algebraic errors"""
    pass class BasisNotSetError(AlgebraError):
    """Raised if the basis or a Hilbert space dimension is unavailable""" class CannotConvertToSLH(AlgebraException):
    """Raised when a circuit algebra object cannot be converted to SLH""" class CannotEliminateAutomatically(AlgebraError):
    """Raised when attempted automatic adiabatic elimination fails.""" class CannotVisualize(AlgebraException):
    """Raised when a circuit cannot be visually represented.""" class IncompatibleBlockStructures(AlgebraError):
    """Raised for invalid block-decomposition

    This is raised when a circuit decomposition into a block-structure is
    requested that is icompatible with the actual block structure of the
    circuit expression.""" class WrongCDimError(AlgebraError):
    """Raised for mismatched channel number in circuit series""" class LocalSpace(HilbertSpace, Expression):
    """Hilbert space for a single degree of freedom.

    Args:
        label (str or int or StrLabel): label (subscript) of the
            Hilbert space
        basis (tuple or None): Set an explicit basis for the Hilbert space
            (tuple of labels for the basis states)
        dimension (int or None): Specify the dimension $n$ of the Hilbert
            space.  This implies a basis numbered from 0  to $n-1$.
        local_identifiers (dict): Mapping of class names of
            :class:`.LocalOperator` subclasses to identifier names. Used e.g.
            'b' instead of the default 'a' for the anihilation operator. This
            can be a dict or a dict-compatible structure, e.g. a list/tuple of
            key-value tuples.
        order_index (int or None): An optional key that determines the
            preferred order of Hilbert spaces. This also changes the order of
            e.g. sums or products of Operators. Hilbert spaces will be ordered
            from left to right be increasing `order_index`; Hilbert spaces
            without an explicit `order_index` are sorted by their label

    A :class:`LocalSpace` fundamentally has a Fock-space like structure,
    in that its basis states may be understood as an "excitation".
    The spectrum can be infinite, with levels labeled by integers 0, 1, ...::

        >>> hs = LocalSpace(label=0)

    or truncated to a finite dimension::

        >>> hs = LocalSpace(0, dimension=5)
        >>> hs.basis_labels
        ('0', '1', '2', '3', '4')

    For finite-dimensional (truncated) Hilbert spaces, we also allow an
    arbitrary alternative labeling of the canonical basis::

        >>> hs = LocalSpace('rydberg', dimension=3, basis=('g', 'e', 'r'))

    """
    _rx_label = re.compile('^[A-Za-z0-9.+-]+(_[A-Za-z0-9().+-]+)?$')

    _basis_label_types = (int, str, FockIndex, FockLabel)
    # acceptable types for labels

    def __init__(
            self, label, *, basis=None, dimension=None, local_identifiers=None,
            order_index=None):

        default_args = []
        if basis is None:
            default_args.append('basis')
        else:
            basis = tuple([str(label) for label in basis])
        if dimension is None:
            default_args.append('dimension')
        else:
            dimension = int(dimension)
        if order_index in [None, float('inf')]:
            default_args.append('order_index')
            order_index = float('inf')  # ensure sort as last
        else:
            order_index = int(order_index)
        if local_identifiers is None:
            default_args.append('local_identifiers')
            local_identifiers = {}
        else:
            local_identifiers = dict(local_identifiers)
        try:
            # we want to normalize the local_identifiers to an arbitrary stable
            # order
            self._sorted_local_identifiers = tuple(
                sorted(tuple(local_identifiers.items())))
        except TypeError:
            # this will happen e.g. if the keys in local_identifier are types
            # instead of class names
            raise TypeError(
                "local_identifier must map class names to identifier strings")

        if not isinstance(label, StrLabel):
            label = str(label)
            if not self._rx_label.match(label):
                raise ValueError(
                    "label '%s' does not match pattern '%s'"
                    % (label, self._rx_label.pattern))
        if basis is None:
            if dimension is not None:
                basis = tuple([str(i) for i in range(dimension)])
        else:
            if dimension is None:
                dimension = len(basis)
            else:
                if dimension != len(basis):
                    raise ValueError("basis and dimension are incompatible")

        self._label = label
        self._order_key = KeyTuple((
            order_index, label, str(dimension), basis,
            self._sorted_local_identifiers))
        self._basis = basis
        self._dimension = dimension
        self._local_identifiers = local_identifiers
        self._order_index = order_index
        self._kwargs = OrderedDict([
            ('basis', self._basis), ('dimension', self._dimension),
            ('local_identifiers', self._sorted_local_identifiers),
            ('order_index', self._order_index)])
        self._minimal_kwargs = self._kwargs.copy()
        for key in default_args:
            del self._minimal_kwargs[key]

        super().__init__(
            label, basis=basis, dimension=dimension,
            local_identifiers=self._sorted_local_identifiers,
            order_index=order_index)

    @classmethod
    def _check_basis_label_type(cls, label_or_index):
        """Every object (BasisKet, LocalSigma) that contains a label or index
        for an eigenstate of some LocalSpace should call this routine to check
        the type of that label or index (or, use
        :meth:`_unpack_basis_label_or_index`"""
        if not isinstance(label_or_index, cls._basis_label_types):
            raise TypeError(
                "label_or_index must be an instance of one of %s; not %s" % (
                    ", ".join([t.__name__ for t in cls._basis_label_types]),
                    label_or_index.__class__.__name__))

    def _unpack_basis_label_or_index(self, label_or_index):
        """return tuple (label, ind) from `label_or_index`

        If `label_or_int` is a :class:`.SymbolicLabelBase` sub-instance, it
        will be stored in the `label` attribute, and the `ind` attribute will
        return the value of the label's :attr:`.FockIndex.fock_index`
        attribute.  No checks are performed for symbolic labels.

        :meth:`_check_basis_label_type` is called on `label_or_index`.

        Raises:
            ValueError: if `label_or_index` is a :class:`str` referencing an
                invalid basis state; or, if `label_or_index` is an :class:`int`
                < 0 or >= the dimension of the Hilbert space
            BasisNotSetError: if `label_or_index` is a :class:`str`, but the
                Hilbert space has no defined basis
            TypeError: if `label_or_int` is not a :class:`str`, :class:`int`,
                or :class:`.SymbolicLabelBase`, or more generally whatever
                types are allowed through the `_basis_label_types` attribute of
                the Hilbert space.
        """
        self._check_basis_label_type(label_or_index)
        if isinstance(label_or_index, str):
            label = label_or_index
            try:
                ind = self.basis_labels.index(label)
                # the above line may also raise BasisNotSetError, which we
                # don't want to catch here
            except ValueError:
                # a less confusing error message:
                raise ValueError(
                    "%r is not one of the basis labels %r"
                    % (label, self.basis_labels))
        elif isinstance(label_or_index, int):
            ind = label_or_index
            if ind < 0:
                raise ValueError("Index %d must be >= 0" % ind)
            if self.has_basis:
                if ind >= self.dimension:
                    raise ValueError(
                        "Index %s must be < the dimension %d of Hilbert "
                        "space %s" % (ind, self.dimension, self))
                label = self.basis_labels[label_or_index]
            else:
                label = str(label_or_index)
        elif isinstance(label_or_index, SymbolicLabelBase):
            label = label_or_index
            try:
                ind = label_or_index.fock_index
            except AttributeError:
                raise TypeError(
                    "label_or_index must define a fock_index attribute in "
                    "order to be used for identifying a level in a Hilbert "
                    "space")
        else:
            raise TypeError(
                "label_or_index must be an int or str, or SymbolicLabelBase, "
                "not %s" % type(label_or_index))
        return label, ind

    @property
    def args(self):
        """List of arguments, consisting only of `label`"""
        return (self._label, )

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

    @property
    def has_basis(self):
        """True if the Hilbert space has a basis"""
        return self._basis is not None

    @property
    def basis_states(self):
        """Yield an iterator over the states (:class:`.BasisKet` instances)
        that form the canonical basis of the Hilbert space

        Raises:
            .BasisNotSetError: if the Hilbert space has no defined basis
        """
        from qnet.algebra.core.state_algebra import BasisKet  # avoid circ. import
        for label in self.basis_labels:
            yield BasisKet(label, hs=self)

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

        Raises:
            .BasisNotSetError: if the Hilbert space has no defined basis
            IndexError: if there is no basis state with the given index
            KeyError: if there is not basis state with the given label
        """
        from qnet.algebra.core.state_algebra import BasisKet  # avoid circ. import
        try:
            return BasisKet(index_or_label, hs=self)
        except ValueError as exc_info:
            if isinstance(index_or_label, int):
                raise IndexError(str(exc_info))
            else:
                raise KeyError(str(exc_info))

    @property
    def basis_labels(self):
        """Tuple of basis labels (strings).

        Raises:
            .BasisNotSetError: if the Hilbert space has no defined basis
        """
        if self._basis is None:
            raise BasisNotSetError(
                "Hilbert space %s has no defined basis" % str(self))
        return self._basis

    @property
    def dimension(self):
        """Dimension of the Hilbert space.

        Raises:
            .BasisNotSetError: if the Hilbert space has no defined basis
        """
        if self._dimension is not None:
            return self._dimension
        else:
            raise BasisNotSetError(
                "Hilbert space %s has no defined basis" % str(self))

    @property
    def kwargs(self):
        return self._kwargs

    @property
    def minimal_kwargs(self):
        return self._minimal_kwargs

    def remove(self, other):
        if other == self:
            return TrivialSpace
        return self

    def intersect(self, other):
        if self in other.local_factors:
            return self
        return TrivialSpace

    @property
    def local_factors(self):
        return (self, )

    def is_strict_subfactor_of(self, other):
        if isinstance(other, ProductSpace) and self in other.operands:
            assert len(other.operands) > 1
            return True
        if other is FullSpace:
            return True
        return False

    def next_basis_label_or_index(self, label_or_index, n=1):
        """Given the label or index of a basis state, return the label/index of
        the next basis state.

        More generally, if `n` is given, return the `n`'th next basis state
        label/index; `n` may also be negative to obtain previous basis state
        labels/indices.

        The return type is the same as the type of `label_or_index`.

        Args:
            label_or_index (int or str or SymbolicLabelBase): If `int`, the
                index of a basis state; if `str`, the label of a basis state
            n (int): The increment

        Raises:
            IndexError: If going beyond the last or first basis state
            ValueError: If `label` is not a label for any basis state in the
                Hilbert space
            .BasisNotSetError: If the Hilbert space has no defined basis
            TypeError: if `label_or_index` is neither a :class:`str` nor an
                :class:`int`, nor a :class:`SymbolicLabelBase`
        """
        if isinstance(label_or_index, int):
            new_index = label_or_index + n
            if new_index < 0:
                raise IndexError("index %d < 0" % new_index)
            if self.has_basis:
                if new_index >= self.dimension:
                    raise IndexError("index %d out of range for basis %s"
                                     % (new_index, self._basis))
            return new_index
        elif isinstance(label_or_index, str):
            label_index = self.basis_labels.index(label_or_index)
            new_index = label_index + n
            if (new_index < 0) or (new_index >= len(self._basis)):
                raise IndexError("index %d out of range for basis %s"
                                 % (new_index, self._basis))
            return self._basis[new_index]
        elif isinstance(label_or_index, SymbolicLabelBase):
            return label_or_index.__class__(expr=label_or_index.expr + n)
        else:
            raise TypeError(
                "Invalid type for label_or_index: %s"
                % label_or_index.__class__.__name__) class ProductSpace(HilbertSpace, Operation):
    """Tensor product of local Hilbert spaces

    >>> hs1 = LocalSpace('1', basis=(0,1))
    >>> hs2 = LocalSpace('2', basis=(0,1))
    >>> hs = hs1 * hs2
    >>> hs.basis_labels
    ('0,0', '0,1', '1,0', '1,1')
    """

    _neutral_element = TrivialSpace
    simplifications = [empty_trivial, assoc, convert_to_spaces, idem,
                       filter_neutral]

    def __init__(self, *local_spaces):
        if len(set(local_spaces)) != len(local_spaces):
            raise ValueError("Nondistinct spaces: %s" % repr(local_spaces))
        try:
            self._dimension = functools.reduce(
                    operator.mul,
                    [ls.dimension for ls in local_spaces], 1)
        except BasisNotSetError:
            self._dimension = None
        # determine the basis labels automatically
        try:
            ls_bases = [ls.basis_labels for ls in local_spaces]
            basis = []
            for label_tuple in cartesian_product(*ls_bases):
                basis.append(",".join([str(l) for l in label_tuple]))
            self._basis = tuple(basis)
        except BasisNotSetError:
            self._basis = None
        op_keys = [space._order_key for space in local_spaces]
        self._order_key = KeyTuple([v for op_key in op_keys for v in op_key])
        super().__init__(*local_spaces)  # Operation __init__

    @classmethod
    def create(cls, *local_spaces):
        if any(local_space is FullSpace for local_space in local_spaces):
            return FullSpace
        return super().create(*local_spaces)

    @property
    def has_basis(self):
        """True if the all the local factors of the `ProductSpace` have a
        defined basis"""
        return self._basis is not None

    @property
    def basis_states(self):
        """Yield an iterator over the states (:class:`.TensorKet` instances)
        that form the canonical basis of the Hilbert space

        Raises:
            .BasisNotSetError: if the Hilbert space has no defined basis
        """
        from qnet.algebra.core.state_algebra import BasisKet, TensorKet
        # importing locally avoids circular import
        ls_bases = [ls.basis_labels for ls in self.local_factors]
        for label_tuple in cartesian_product(*ls_bases):
            yield TensorKet(
                *[BasisKet(label, hs=ls) for (ls, label)
                  in zip(self.local_factors, label_tuple)])

    @property
    def basis_labels(self):
        """Tuple of basis labels. Each basis label consists of the labels of
        the :class:`.BasisKet` states that factor the basis state, separated by
        commas.

        Raises:
            .BasisNotSetError: if the Hilbert space has no defined basis
        """
        if self._basis is None:
            raise BasisNotSetError(
                "Hilbert space %s has no defined basis" % str(self))
        return self._basis

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

        Raises:
            .BasisNotSetError: if the Hilbert space has no defined basis
            IndexError: if there is no basis state with the given index
            KeyError: if there is not basis state with the given label
        """
        from qnet.algebra.core.state_algebra import BasisKet, TensorKet
        if isinstance(index_or_label, int):  # index
            ls_bases = [ls.basis_labels for ls in self.local_factors]
            label_tuple = list(cartesian_product(*ls_bases))[index_or_label]
            try:
                return TensorKet(
                    *[BasisKet(label, hs=ls) for (ls, label)
                        in zip(self.local_factors, label_tuple)])
            except ValueError as exc_info:
                raise IndexError(str(exc_info))
        else:  # label
            local_labels = index_or_label.split(",")
            if len(local_labels) != len(self.local_factors):
                raise KeyError(
                    "label %s for Hilbert space %s must be comma-separated "
                    "concatenation of local labels" % (index_or_label, self))
            try:
                return TensorKet(
                    *[BasisKet(label, hs=ls) for (ls, label)
                      in zip(self.local_factors, local_labels)])
            except ValueError as exc_info:
                raise KeyError(str(exc_info))

    @property
    def dimension(self):
        """Dimension of the Hilbert space.

        Raises:
            .BasisNotSetError: if the Hilbert space has no defined basis
        """
        if self._dimension is not None:
            return self._dimension
        else:
            raise BasisNotSetError(
                "Hilbert space %s has no defined basis" % str(self))

    def remove(self, other):
        """Remove a particular factor from a tensor product space."""
        if other is FullSpace:
            return TrivialSpace
        if other is TrivialSpace:
            return self
        if isinstance(other, ProductSpace):
            oops = set(other.operands)
        else:
            oops = {other}
        return ProductSpace.create(
                *sorted(set(self.operands).difference(oops)))

    @property
    def local_factors(self):
        """The :class:`LocalSpace` instances that make up the product"""
        return self.operands

    @classmethod
    def order_key(cls, obj):
        """Key by which operands are sorted"""
        assert isinstance(obj, HilbertSpace)
        return obj._order_key

    def intersect(self, other):
        """Find the mutual tensor factors of two Hilbert spaces."""
        if other is FullSpace:
            return self
        if other is TrivialSpace:
            return TrivialSpace
        if isinstance(other, ProductSpace):
            other_ops = set(other.operands)
        else:
            other_ops = {other}
        return ProductSpace.create(
                *sorted(set(self.operands).intersection(other_ops)))

    def is_strict_subfactor_of(self, other):
        """Test if a space is included within a larger tensor product space.
        Not ``True`` if ``self == other``."""
        if isinstance(other, ProductSpace):
            return set(self.operands) < set(other.operands)
        if other is FullSpace:
            return True
        return False class Matrix(Expression):
    """Matrix of Expressions

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

                (M.conjugate() + M) / 2

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

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

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

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

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

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

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

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

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

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

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

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

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

    dag = adjoint

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        return self.element_wise(element_simplify) def block_matrix(A, B, C, D):
    r"""Generate the operator matrix with quadrants

    .. math::

       \begin{pmatrix} A B \\ C D \end{pmatrix}

    Args:
        A (Matrix): Matrix of shape ``(n, m)``
        B (Matrix): Matrix of shape ``(n, k)``
        C (Matrix): Matrix of shape ``(l, m)``
        D (Matrix): Matrix of shape ``(l, k)``

    Returns:
        Matrix: The combined block matrix ``[[A, B], [C, D]]``.
    """
    return vstackm((hstackm((A, B)), hstackm((C, D)))) def identity_matrix(N):
    """Generate the N-dimensional identity matrix.

    Args:
        N (int): Dimension

    Returns:
        Matrix: Identity matrix in N dimensions

    """
    return diagm(np_ones(N, dtype=int)) def permutation_matrix(permutation):
    r"""Return orthogonal permutation matrix for permutation tuple

    Return an orthogonal permutation matrix :math:`M_\sigma`
    for a permutation :math:`\sigma` defined by the image tuple
    :math:`(\sigma(1), \sigma(2),\dots \sigma(n))`,
    such that

    .. math::

        M_\sigma \vec{e}_i = \vec{e}_{\sigma(i)}

    where :math:`\vec{e}_k` is the k-th standard basis vector.
    This definition ensures a composition law:

    .. math::

        M_{\sigma \cdot \tau} = M_\sigma M_\tau.

    The column form of :math:`M_\sigma` is thus given by

    .. math::

        M = (
            \vec{e}_{\sigma(1)},
            \vec{e}_{\sigma(2)},
            \dots \vec{e}_{\sigma(n)}).

    Args:
        permutation (tuple): A permutation image tuple (zero-based indices!)
    """
    assert check_permutation(permutation)
    n = len(permutation)
    op_matrix = np_zeros((n, n), dtype=int)
    for i, j in enumerate(permutation):
        op_matrix[j, i] = 1
    return Matrix(op_matrix) def vstackm(matrices):
    """Generalizes `numpy.vstack` to :class:`Matrix` objects."""
    arr = np_vstack(tuple(m.matrix for m in matrices))
    #    print(tuple(m.matrix.dtype for m in matrices))
    #    print(arr.dtype)
    return Matrix(arr) def zerosm(shape, *args, **kwargs):
    """Generalizes ``numpy.zeros`` to :py:class:`Matrix` objects."""
    return Matrix(np_zeros(shape, *args, **kwargs)) class IdentityOperator(Operator, metaclass=Singleton):
    """``IdentityOperator`` constant (singleton) object."""

    _order_index = 2

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

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

    def _diff(self, sym):
        return ZeroOperator

    def _adjoint(self):
        return self

    def _pseudo_inverse(self):
        return self def LocalProjector(j, *, hs):
    """A projector onto a specific level of a :class:`.LocalSpace`

    Args:
        j (int or str): The label or index identifying the state onto which
            is projected
        hs (HilbertSpace): The Hilbert space on which the operator acts
    """
    return LocalSigma(j, j, hs=hs) class LocalSigma(LocalOperator):
    r'''Level flip operator between two levels of a :class:`.LocalSpace`

    .. math::

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        .. math::

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

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

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

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

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

    def _pseudo_inverse(self):
        return self._adjoint() class 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 OperatorSymbol(QuantumSymbol, Operator):
    """Symbolic operator

    See :class:`.QuantumSymbol`.
    """
    def _pseudo_inverse(self):
        return PseudoInverse(self) class ScalarTimesOperator(Operator, ScalarTimesQuantumExpression):
    """Product of a :class:`.Scalar` coefficient and an :class:`Operator`"""

    _rules = OrderedDict()
    simplifications = [match_replace, ]

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

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

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

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

    def _adjoint(self):
        return ScalarTimesOperator(self.coeff.conjugate(), self.term.adjoint()) class ZeroOperator(Operator, metaclass=Singleton):
    """``ZeroOperator`` constant (singleton) object."""

    _order_index = 2

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

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

    def _diff(self, sym):
        return self

    def _adjoint(self):
        return self

    def _pseudo_inverse(self):
        return self def adjoint(obj):
    """Return the adjoint of an obj."""
    try:
        return obj.adjoint()
    except AttributeError:
        return obj.conjugate() def get_coeffs(expr, expand=False, epsilon=0.):
    """Create a dictionary with all Operator terms of the expression
    (understood as a sum) as keys and their coefficients as values.

    The returned object is a defaultdict that return 0. if a term/key
    doesn't exist.

    Args:
        expr: The operator expression to get all coefficients from.
        expand: Whether to expand the expression distributively.
        epsilon: If non-zero, drop all Operators with coefficients that have
            absolute value less than epsilon.

    Returns:
        dict: A dictionary ``{op1: coeff1, op2: coeff2, ...}``
    """
    if expand:
        expr = expr.expand()
    ret = defaultdict(int)
    operands = expr.operands if isinstance(expr, OperatorPlus) else [expr]
    for e in operands:
        c, t = _coeff_term(e)
        try:
            if abs(complex(c)) < epsilon:
                continue
        except TypeError:
            pass
        ret[t] += c
    return ret class BadPermutationError(ValueError):
    """Can be raised to signal that a permutation does not pass the
    :py:func:check_permutation test."""
    pass def block_perm_and_perms_within_blocks(permutation, block_structure):
    """Decompose a permutation into a block permutation and into permutations
    acting within each block.

    :param permutation: The overall permutation to be factored.
    :type permutation: tuple
    :param block_structure: The channel dimensions of the blocks
    :type block_structure: tuple
    :return: ``(block_permutation, permutations_within_blocks)``
     Where ``block_permutations`` is an image tuple for a permutation of the block indices
     and ``permutations_within_blocks`` is a list of image tuples for the permutations of the channels
     within each block
    :rtype: tuple
    """
    nblocks = len(block_structure)

    offsets = [sum(block_structure[:k]) for k in range(nblocks)]
    images = [permutation[offset: offset + length] for (offset, length) in zip(offsets, block_structure)]

    images_mins = list(map(min, images))


    key_block_perm_inv = lambda block_index: images_mins[block_index]

    block_perm_inv = tuple(sorted(range(nblocks), key = key_block_perm_inv))
    # print(images_mins)
    # print(permutation, block_structure, "-->", block_perm, invert_permutation(block_perm))
    block_perm = invert_permutation(block_perm_inv)

    assert images_mins[block_perm_inv[0]] == min(images_mins)
    assert images_mins[block_perm_inv[-1]] == max(images_mins)

    # block_perm = tuple(invert_permutation(block_perm_inv))

    perms_within_blocks = []
    for (offset, length, image) in zip(offsets, block_structure, images):
        block_key = lambda elt_index: image[elt_index]
        within_inv = sorted(range(length), key = block_key)
        within = invert_permutation(tuple(within_inv))
        assert permutation[within_inv[0] + offset] == min(image)
        assert permutation[within_inv[-1] + offset] == max(image)
        perms_within_blocks.append(within)

    return block_perm, perms_within_blocks def check_permutation(permutation):
    """Verify that a tuple of permutation image points ``(sigma(1), sigma(2),
    ..., sigma(n))`` is a valid permutation, i.e. each number from 0 and n-1
    occurs exactly once. I.e. the following **set**-equality must hold:

        ``{sigma(1), sigma(2), ..., sigma(n)} == {0, 1, 2, ... n-1}``

    :param permutation: Tuple of permutation image points
    :type permutation: tuple
    :rtype: bool
    """
    return list(sorted(permutation)) == list(range(len(permutation))) def full_block_perm(block_permutation, block_structure):
    """
    Extend a permutation of blocks to a permutation for the internal signals of all blocks.
    E.g., say we have two blocks of sizes ('block structure') ``(2, 3)``,
    then a block permutation that switches the blocks would be given by the image tuple ``(1,0)``.
    However, to get a permutation of all 2+3 = 5 channels that realizes that block permutation we would need
    ``(2, 3, 4, 0, 1)``

    :param block_permutation: permutation image tuple of block indices
    :type block_permutation: tuple
    :param block_structure: The block channel dimensions, block structure
    :type block_structure: tuple
    :return: A single permutation for all channels of all blocks.
    :rtype: tuple
    """
    fblockp = []
    bp_inv = invert_permutation(block_permutation)
    for k, block_length in enumerate(block_structure):
        p_k = block_permutation[k]
        offset = sum([block_structure[bp_inv[j]] for j in range(p_k)])
        fblockp += range(offset, offset + block_length)

    assert sorted(fblockp) == list(range(sum(block_structure)))

    return tuple(fblockp) def invert_permutation(permutation):
    """Compute the image tuple of the inverse permutation.

    :param permutation: A valid (cf. :py:func:check_permutation) permutation.
    :return: The inverse permutation tuple
    :rtype: tuple
    """
    return tuple([permutation.index(p) for p in range(len(permutation))]) def permutation_to_block_permutations(permutation):
    """If possible, decompose a permutation into a sequence of permutations
    each acting on individual ranges of the full range of indices.
    E.g.

        ``(1,2,0,3,5,4) --> (1,2,0) [+] (0,2,1)``

    :param permutation: A valid permutation image tuple ``s = (s_0,...s_n)`` with ``n > 0``
    :type permutation: tuple
    :return: A list of permutation tuples ``[t = (t_0,...,t_n1), u = (u_0,...,u_n2),..., z = (z_0,...,z_nm)]`` such that ``s = t [+] u [+] ... [+] z``
    :rtype: list of tuples
    :raise: ValueError
    """
    if len(permutation) == 0 or not check_permutation(permutation):
        raise BadPermutationError()

    cycles = permutation_to_disjoint_cycles(permutation)

    if len(cycles) == 1:
        return (permutation,)
    current_block_start = cycles[0][0]
    current_block_end = max(cycles[0])
    current_block_cycles = [cycles[0]]
    res_permutations = []
    for c in cycles[1:]:
        if c[0] > current_block_end:
            res_permutations.append(permutation_from_disjoint_cycles(current_block_cycles, current_block_start))
            assert sum(map(len, current_block_cycles)) == current_block_end - current_block_start + 1
            current_block_start = c[0]
            current_block_end = max(c)
            current_block_cycles = [c]
        else:
            current_block_cycles.append(c)
            if max(c) > current_block_end:
                current_block_end = max(c)

    res_permutations.append(permutation_from_disjoint_cycles(current_block_cycles, current_block_start))

    assert sum(map(len, current_block_cycles)) == current_block_end - current_block_start + 1
    assert sum(map(len, res_permutations)) == len(permutation)

    return res_permutations 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 os
import re
import numpy as np
import sympy
        import qnet.visualization.circuit_pyx as circuit_visualization
from abc import ABCMeta, abstractmethod
from collections import OrderedDict
from functools import reduce
from sympy import I
from sympy import Matrix as SympyMatrix
from sympy import symbols, sympify
from .abstract_algebra import (
    Expression, Operation, substitute, )
from .algebraic_properties import (
    assoc, check_cdims, filter_neutral, filter_cid, match_replace,
    match_replace_binary)
from .exceptions import (
    AlgebraError, BasisNotSetError, CannotConvertToSLH,
    CannotEliminateAutomatically, CannotVisualize, IncompatibleBlockStructures,
    WrongCDimError)
from .hilbert_space_algebra import LocalSpace, ProductSpace
from .matrix_algebra import (
    Matrix, block_matrix, identity_matrix, permutation_matrix,
    vstackm, zerosm)
from .operator_algebra import (
    IdentityOperator, LocalProjector, LocalSigma, Operator,
    OperatorPlus, OperatorSymbol, ScalarTimesOperator, ZeroOperator, adjoint,
    get_coeffs)
from ...utils.permutations import (
    BadPermutationError, block_perm_and_perms_within_blocks, check_permutation,
    full_block_perm, invert_permutation, permutation_to_block_permutations, )
from ...utils.singleton import Singleton, singleton_object
        from IPython.display import Image, display
        from tempfile import gettempdir
        from time import time, sleep
        from qnet.algebra.library.circuit_components import (
                CoherentDriveCC as Displace_cc)
        from sympy.core.numbers import ComplexInfinity, Infinity
        from qnet.algebra.core.super_operator_algebra import liouvillian
    from qnet.algebra.library.fock_operators import Create, Destroy
            n           := self.cdim
            in_im       := self.permutation[in_port]
            m_{k->l}    := map_signals_circuit({k:l}, n)

        solve the equation (I) containing ``self``::

            self << m_{(n-1) -> in_port}
                == m_{(n-1) -> in_im} << (red_self + cid(1))          (I)

        for the (n-1) channel CPermutation ``red_self``.
        Return in_im, red_self.

        This is useful when ``self`` is the RHS in a SeriesProduct Object that
        is within a Feedback loop as it allows to extract the feedback channel
        from the permutation and moving the remaining part of the permutation
        (``red_self``) outside of the feedback loop.

        :param int in_port: The index for which to factor.
        """
        n = self.cdim
        if not (0 <= in_port < n):
            raise Exception
        in_im = self.permutation[in_port]
        # (I) is equivalent to
        #       m_{in_im -> (n-1)} <<  self << m_{(n-1) -> in_port}
        #           == (red_self + cid(1))     (I')
        red_self_plus_cid1 = (map_channels({in_im: (n - 1)}, n) <<
                              self <<
                              map_channels({(n - 1): in_port}, n))
        if isinstance(red_self_plus_cid1, CPermutation):

            #make sure we can factor
            assert red_self_plus_cid1.permutation[(n - 1)] == (n - 1)

            #form reduced permutation object
            red_self = CPermutation.create(red_self_plus_cid1.permutation[:-1])

            return in_im, red_self
        else:
            # 'red_self_plus_cid1' must be the identity for n channels.
            # Actually, this case can only occur
            # when self == m_{in_port ->  in_im}

            return in_im, circuit_identity(n - 1)

    def _factor_lhs(self, out_port):
        """With::

            n           := self.cdim
            out_inv     := invert_permutation(self.permutation)[out_port]
            m_{k->l}    := map_signals_circuit({k:l}, n)

        solve the equation (I) containing ``self``::

            m_{out_port -> (n-1)} << self
                == (red_self + cid(1)) << m_{out_inv -> (n-1)}           (I)

        for the (n-1) channel CPermutation ``red_self``.
        Return out_inv, red_self.

        This is useful when 'self' is the LHS in a SeriesProduct Object that is
        within a Feedback loop as it allows to extract the feedback channel
        from the permutation and moving the remaining part of the permutation
        (``red_self``) outside of the feedback loop.
        """
        n = self.cdim
        assert (0 <= out_port < n)
        out_inv = self.permutation.index(out_port)

        # (I) is equivalent to
        #       m_{out_port -> (n-1)} <<  self << m_{(n-1)
        #           -> out_inv} == (red_self + cid(1))     (I')

        red_self_plus_cid1 = (map_channels({out_port: (n - 1)}, n) <<
                              self <<
                              map_channels({(n - 1): out_inv}, n))

        if isinstance(red_self_plus_cid1, CPermutation):

            #make sure we can factor
            assert red_self_plus_cid1.permutation[(n - 1)] == (n - 1)

            #form reduced permutation object
            red_self = CPermutation.create(red_self_plus_cid1.permutation[:-1])

            return out_inv, red_self
        else:
            # 'red_self_plus_cid1' must be the identity for n channels.
            # Actually, this case can only occur
            # when self == m_{in_port ->  in_im}

            return out_inv, circuit_identity(n - 1)


@singleton_object
class CIdentity(Circuit, Expression, metaclass=Singleton):
    """Single pass-through channel; neutral element of :class:`SeriesProduct`
    """

    _cdim = 1

    @property
    def cdim(self):
        """Dimension of circuit"""
        return self._cdim

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

    def _toSLH(self):
        return SLH(Matrix([[1]]), Matrix([[0]]), 0)

    def _creduce(self):
        return self

    def _series_inverse(self):
        return self



@singleton_object