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 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 TrivialSpace(HilbertSpace, Expression, metaclass=Singleton):
    """The 'nullspace', i.e. a one dimensional Hilbert space, which is a factor
    space of every other Hilbert space.

    This is the Hilbert space of scalars.
    """

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

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

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

    @property
    def dimension(self):
        return 1

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

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

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

        All other indices or labels raise an exception.

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

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

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

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

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

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

    @property
    def label(self):
        return "null" class Operation(Expression, metaclass=ABCMeta):
    """Base class for "operations"

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

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

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

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

    @property
    def args(self):
        """Alias for operands"""
        return self._operands class 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 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 derivative_via_diff(cls, ops, kwargs):
    """Implementation of the :meth:`QuantumDerivative.create` interface via the
    use of :meth:`QuantumExpression._diff`.

    Thus, by having :meth:`.QuantumExpression.diff` delegate to
    :meth:`.QuantumDerivative.create`, instead of
    :meth:`.QuantumExpression._diff` directly, we get automatic caching of
    derivatives
    """
    assert len(ops) == 1
    op = ops[0]
    derivs = kwargs['derivs']
    vals = kwargs['vals']
    # both `derivs` and `vals` are guaranteed to be tuples, via the conversion
    # that's happening in `QuantumDerivative.create`
    for (sym, n) in derivs:
        if sym.free_symbols.issubset(op.free_symbols):
            for k in range(n):
                op = op._diff(sym)
        else:
            return op.__class__._zero
    if vals is not None:
        try:
            # for QuantumDerivative instance
            return op.evaluate_at(vals)
        except AttributeError:
            # for explicit Expression
            return op.substitute(vals)
    else:
        return op class IndexedSum(Operation, metaclass=ABCMeta):
    """Base class for indexed sums"""

    _expanded_cls = None  # must be set by subclasses
    _expand_limit = 1000

    def __init__(self, term, *ranges):
        self._term = term
        self.ranges = tuple(ranges)
        for r in self.ranges:
            if not isinstance(r, IndexRangeBase):
                # We need this type check to that we can use attr.astuple below
                raise TypeError(
                    "Every range must be an instance of IndexRangeBase")

        index_symbols = set([r.index_symbol for r in ranges])
        if len(index_symbols) != len(self.ranges):
            raise ValueError(
                "ranges %s must have distinct index_symbols" % repr(ranges))
        super().__init__(term, ranges=ranges)

    @property
    def term(self):
        return self._term

    @property
    def operands(self):
        return (self._term, )

    @property
    def args(self):
        return tuple([self._term, *self.ranges])

    @property
    def variables(self):
        """List of the dummy (index) variable symbols

        See also :property:`bound_symbols` for a set of the same symbols
        """
        return [r.index_symbol for r in self.ranges]

    @property
    def bound_symbols(self):
        """Set of bound variables, i.e. the index variable symbols

        See also :property:`variables` for an ordered list of the same symbols
        """
        return set(self.variables)

    @property
    def free_symbols(self):
        """Set of all free symbols"""
        return set([
            sym for sym in self.term.free_symbols
            if sym not in self.bound_symbols])

    @property
    def kwargs(self):
        return {}

    @property
    def terms(self):
        """Iterator over the terms of the sum

        Yield from the (possibly) infinite list of terms of the indexed sum, if
        the sum was written out explicitly. Each yielded term in an instance of
        :class:`.Expression`
        """
        from qnet.algebra.core.scalar_algebra import ScalarValue
        for mapping in yield_from_ranges(self.ranges):
            term = self.term.substitute(mapping)
            if isinstance(term, ScalarValue._val_types):
                term = ScalarValue.create(term)
            assert isinstance(term, Expression)
            yield term

    def __len__(self):
        length = 1
        for ind_range in self.ranges:
            try:
                length *= len(ind_range)
            except TypeError:
                raise InfiniteSumError(
                    "Cannot determine length from non-finite ranges")
        return length

    def doit(
            self, classes=None, recursive=True, indices=None, max_terms=None,
            **kwargs):
        """Write out the indexed sum explicitly

        If `classes` is None or :class:`IndexedSum` is in `classes`,
        (partially) write out the indexed sum in to an explicit sum of terms.
        If `recursive` is True, write out each of the new sum's summands by
        calling its :meth:`doit` method.

        Args:
            classes (None or list): see :meth:`.Expression.doit`
            recursive (bool): see :meth:`.Expression.doit`
            indices (list): List of :class:`IdxSym` indices for which the sum
                should be expanded. If `indices` is a subset of the indices
                over which the sum runs, it will be partially expanded. If not
                given, expand the sum completely
            max_terms (int): Number of terms after which to truncate the sum.
                This is particularly useful for infinite sums. If not given,
                expand all terms of the sum. Cannot be combined with `indices`
            kwargs: keyword arguments for recursive calls to
                :meth:`doit`. See :meth:`.Expression.doit`
        """
        return super().doit(
            classes, recursive, indices=indices, max_terms=max_terms, **kwargs)

    def _doit(self, **kwargs):
        indices = kwargs.get('indices', None)
        max_terms = kwargs.get('max_terms', None)
        if indices is None:
            return self._doit_full(max_terms=max_terms)
        else:
            if max_terms is not None:
                raise ValueError(
                    "max_terms is incompatible with summing over specific "
                    "indices")
            return self._doit_over_indices(indices)

    def _doit_full(self, max_terms=None):
        res = None
        if max_terms is None:
            len(self)  # side-effect: raise InfiniteSumError
        else:
            if max_terms > self._expand_limit:
                raise ValueError(
                    "max_terms = %s must be smaller than the limit %s"
                    % (max_terms, self._expand_limit))
        for i, term in enumerate(self.terms):
            if max_terms is not None:
                if i >= max_terms:
                    break
            if res is None:
                res = term
            else:
                res += term
            if i > self._expand_limit:
                raise InfiniteSumError(
                    "Cannot expand %s: more than %s terms"
                    % (self, self._expand_limit))
        return res

    def _doit_over_indices(self, indices):
        if len(indices) == 0:
            return self
        ind_sym, *indices = indices
        if not isinstance(ind_sym, IdxSym):
            ind_sym = IdxSym(ind_sym)
        selected_range = None
        other_ranges = []
        for index_range in self.ranges:
            if index_range.index_symbol == ind_sym:
                selected_range = index_range
            else:
                other_ranges.append(index_range)
        if selected_range is None:
            return self
        res_term = None
        for i, mapping in enumerate(selected_range.iter()):
            res_summand = self.term.substitute(mapping)
            if res_term is None:
                res_term = res_summand
            else:
                res_term += res_summand
            if i > self._expand_limit:
                raise InfiniteSumError(
                    "Cannot expand %s: more than %s terms"
                    % (self, self._expand_limit))
        if len(other_ranges) == 0:
            res = res_term
        else:
            res = self.__class__.create(res_term, *other_ranges)
            res = res._doit_over_indices(indices=indices)
        return res

    def make_disjunct_indices(self, *others):
        """Return a copy with modified indices to ensure disjunct indices with
        `others`.

        Each element in `others` may be an index symbol (:class:`.IdxSym`),
        a index-range object (:class:`.IndexRangeBase`) or list of index-range
        objects, or an indexed operation (something with a ``ranges`` attribute)

        Each index symbol is primed until it does not match any index symbol in
        `others`.
        """
        new = self
        other_index_symbols = set()
        for other in others:
            try:
                if isinstance(other, IdxSym):
                    other_index_symbols.add(other)
                elif isinstance(other, IndexRangeBase):
                    other_index_symbols.add(other.index_symbol)
                elif hasattr(other, 'ranges'):
                    other_index_symbols.update(
                        [r.index_symbol for r in other.ranges])
                else:
                    other_index_symbols.update(
                        [r.index_symbol for r in other])
            except AttributeError:
                raise ValueError(
                    "Each element of other must be an an instance of IdxSym, "
                    "IndexRangeBase, an object with a `ranges` attribute "
                    "with a list of IndexRangeBase instances, or a list of"
                    "IndexRangeBase objects.")
        for r in self.ranges:
            index_symbol = r.index_symbol
            modified = False
            while index_symbol in other_index_symbols:
                modified = True
                index_symbol = index_symbol.incr_primed()
            if modified:
                new = new._substitute(
                    {r.index_symbol: index_symbol}, safe=True)
        return new

    def __mul__(self, other):
        if isinstance(other, IndexedSum):
            other = other.make_disjunct_indices(self)
            new_ranges = self.ranges + other.ranges
            return self.__class__.create(self.term * other.term, *new_ranges)
        try:
            return super().__mul__(other)
        except AttributeError:
            return NotImplemented

    def __rmul__(self, other):
        if isinstance(other, IndexedSum):
            self = self.make_disjunct_indices(other)
            new_ranges = other.ranges + self.ranges
            return self.__class__.create(other.term * self.term, *new_ranges)
        try:
            return super().__rmul__(other)
        except AttributeError:
            return NotImplemented class DisjunctCommutativeHSOrder():
    """Auxiliary class that generates the correct pseudo-order relation for
    operator products.  Only operators acting on disjoint Hilbert spaces
    are commuted to reflect the order the local factors have in the total
    Hilbert space. I.e., ``sorted(factors, key=DisjunctCommutativeHSOrder)``
    achieves this ordering.
    """

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

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

    def __lt__(self, other):
        if self.trivial and other.trivial:
            return self._op_order < other._op_order
        else:
            if self.space.isdisjoint(other.space):
                return self._space_order < other._space_order
        return None  # no ordering class FullCommutativeHSOrder():
    """Auxiliary class that generates the correct pseudo-order relation for
    operator sums.  Operators are first ordered by their Hilbert space, then by
    their order-key; ``sorted(factors, key=FullCommutativeHSOrder)``
    achieves this ordering.
    """

    def __init__(self, op, space_order=None, op_order=None):
        self.op = op
        self.space = op.space
        if space_order is None:
            self._space_order = self.space._order_key
        else:
            self._space_order = space_order
        if op_order is None:
            self._op_order = expr_order_key(op)
        else:
            self._op_order = op_order

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

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

    def __repr__(self):
        return self.__class__.__name__ + tuple.__repr__(self) class SymbolicLabelBase(metaclass=ABCMeta):
    """Base class for symbolic labels

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

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

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

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

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

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

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

    @property
    def free_symbols(self):
        """Free symbols in the expression describing the label"""
        return self.expr.free_symbols class IndexOverList(IndexRangeBase):
    """Index over a list of explicit values

    Args:
        index_symbol (IdxSym): The symbol iterating over the value
        values (list): List of values for the index
    """
    values = attr.ib(convert=tuple)

    def iter(self):
        for val in self.values:
            yield {self.index_symbol: val}

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

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

    def substitute(self, var_map):
        new_index_symbol = var_map.get(self.index_symbol, self.index_symbol)
        new_values = tuple(
            [var_map.get(element, element) for element in self.values])
        return self.__class__(index_symbol=new_index_symbol, values=new_values)

    def piecewise_one(self, expr):
        return Piecewise(
            (1, sympy.FiniteSet(*self.values).contains(expr)),
            (0, True)) class IndexOverFockSpace(IndexRangeBase):
    """Index range over the integer indices of a :class:`.LocalSpace` basis

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

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

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

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

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

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

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

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

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

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

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

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

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

    def piecewise_one(self, expr):
        return Piecewise(
            (1, sympy.FiniteSet(*list(self.range)).contains(expr)),
            (0, True))
import re
import sympy
from abc import ABCMeta, abstractmethod
from collections import defaultdict, OrderedDict
from itertools import product as cartesian_product
from sympy import Symbol, sympify
from .hilbert_space_algebra import ProductSpace, LocalSpace, TrivialSpace
from .abstract_algebra import Operation, Expression, substitute
from .algebraic_properties import derivative_via_diff
from .indexed_operations import IndexedSum
from ...utils.ordering import (
    DisjunctCommutativeHSOrder, FullCommutativeHSOrder, KeyTuple, )
from ...utils.indices import (
    SymbolicLabelBase, IndexOverList, IndexOverFockSpace, IndexOverRange)
        from qnet.algebra.core.scalar_algebra import is_scalar, ScalarValue
        from qnet.algebra.core.scalar_algebra import is_scalar
        from qnet.algebra.core.scalar_algebra import ScalarValue
        from qnet.algebra.core.scalar_algebra import Scalar, ScalarValue
        from qnet.algebra.core.scalar_algebra import Scalar, ScalarValue
        from qnet.printing.asciiprinter import QnetAsciiDefaultPrinter
        from qnet.printing.asciiprinter import QnetAsciiDefaultPrinter
        from qnet.algebra.core.scalar_algebra import is_scalar
        from qnet.algebra.core.scalar_algebra import is_scalar
    from qnet.algebra.core.hilbert_space_algebra import LocalSpace
    from qnet.algebra.core.scalar_algebra import ScalarValue
    from qnet.algebra.library.spin_algebra import SpinSpace
    from qnet.algebra.core.scalar_algebra import Zero
                    in self.operand.substitute(dummy_map).free_symbols
                    if not isinstance(sym, sympy.Dummy)}
                for val in self._vals.values():
                    self._free_symbols.update(val.free_symbols)
        return self._free_symbols

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

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

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

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

    def _adjoint(self):
        return self.__class__(
            self.operand.adjoint(), derivs=self._derivs, vals=self._vals)


class QuantumIndexedSum(IndexedSum, SingleQuantumOperation, metaclass=ABCMeta):
    """Base class for indexed sums"""

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

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

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

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

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

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

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

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

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

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


def _sum_over_list(term, idx, values):
    return IndexOverList(idx, values)


def _sum_over_range(term, idx, start_from, to, step=1):
    return IndexOverRange(idx, start_from=start_from, to=to, step=step)


def _sum_over_fockspace(term, idx, hs=None):
    if hs is None:

        return IndexOverFockSpace(idx, hs=term.space)