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 IdentityOperator(Operator, metaclass=Singleton):
    """``IdentityOperator`` constant (singleton) object."""

    _order_index = 2

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

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

    def _diff(self, sym):
        return ZeroOperator

    def _adjoint(self):
        return self

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

    _order_index = 2

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

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

    def _diff(self, sym):
        return self

    def _adjoint(self):
        return self

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    .. math::

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        .. math::

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

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

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

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

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

    def _pseudo_inverse(self):
        return self._adjoint() class 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 OperatorTimes(QuantumTimes, Operator):
    """Product of operators

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

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

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

    _rules = OrderedDict()
    simplifications = [match_replace, ]

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

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

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

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

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

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

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

    It is defined via the relationship

    .. math::

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

    def _expand(self):
        return self

    def _pseudo_inverse(self):
        return self.operand

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

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

    .. math::

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

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

    def _expand(self):
        return self

    def _adjoint(self):
        return Adjoint(self)

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

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

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

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

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

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

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

    Printers should represent this operator with the default identifier::

        >>> Jz._identifier
        'J_z'

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

    _identifier = 'J_z'

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

    def _adjoint(self):
        return self

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

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

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

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

    Printers should represent this operator with the default identifier::

        >>> Jplus._identifier
        'J_+'

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

    _identifier = 'J_+'

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

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

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

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

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

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

    Printers should represent this operator with the default identifier::

        >>> Jminus._identifier
        'J_-'

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

    _identifier = 'J_-'

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

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

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

    It obeys the bosonic commutation relation::

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

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

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

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

    def _pseudo_inverse(self):
        return PseudoInverse(self)

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

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

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

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

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

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

    def _pseudo_inverse(self):
        return PseudoInverse(self)

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

    .. math::

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

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

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

    Printers should represent this operator with the default identifier::

        >>> Phase._identifier
        'Phase'

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

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

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

    .. math::

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

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

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

    Printers should represent this operator with the default identifier::

        >>> Displace._identifier
        'D'

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

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

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

    .. math::

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

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

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

    Printers should represent this operator with the default identifier::

        >>> Squeeze._identifier
        'Squeeze'

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

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

    _pseudo_inverse = _adjoint
import sympy
from sympy.physics.quantum import TensorProduct as tensor
from ..algebra.core.abstract_algebra import Operation
from ..algebra.core.operator_algebra import (
    IdentityOperator, ZeroOperator, LocalOperator, LocalSigma, Operator,
    OperatorPlus, OperatorTimes, ScalarTimesOperator,
    Adjoint, PseudoInverse, NullSpaceProjector)
from ..algebra.library.spin_algebra import Jz, Jplus, Jminus
from ..algebra.library.fock_operators import (
    Destroy, Create, Phase, Displace, Squeeze)

__all__ = ['convert_to_sympy_matrix']
__private__ = ['SympyCreate', 'basis_state']


def basis_state(i, n):
    """``n x 1`` `sympy.Matrix` representing the `i`'th eigenstate of an
    `n`-dimensional Hilbert space (`i` >= 0)"""
    v = sympy.zeros(n, 1)
    v[i] = 1
    return v


def SympyCreate(n):
    """Creation operator for a Hilbert space of dimension `n`, as an instance
    of `sympy.Matrix`"""
    a = sympy.zeros(n)
    for i in range(1, n):
        a += sympy.sqrt(i) * basis_state(i, n) * basis_state(i-1, n).H
    return a


def convert_to_sympy_matrix(expr, full_space=None):
    """Convert a QNET expression to an explicit ``n x n`` instance of
    `sympy.Matrix`, where ``n`` is the dimension of `full_space`. The entries
    of the matrix may contain symbols.

    Parameters:
        expr: a QNET expression
        full_space (qnet.algebra.hilbert_space_algebra.HilbertSpace): The
            Hilbert space in which `expr` is defined. If not given,
            ``expr.space`` is used. The Hilbert space must have a well-defined
            basis.

    Raises:
        qnet.algebra.hilbert_space_algebra.BasisNotSetError: if `full_space`
            does not have a defined basis
        ValueError: if `expr` is not in `full_space`, or if `expr` cannot be
            converted.
    """
    if full_space is None:
        full_space = expr.space
    if not expr.space.is_tensor_factor_of(full_space):
        raise ValueError("expr must be in full_space")
    if expr is IdentityOperator:
        return sympy.eye(full_space.dimension)
    elif expr is ZeroOperator:
        return 0
    elif isinstance(expr, LocalOperator):
        n = full_space.dimension
        if full_space != expr.space:
            all_spaces = full_space.local_factors
            own_space_index = all_spaces.index(expr.space)
            factors = [sympy.eye(s.dimension)
                       for s in all_spaces[:own_space_index]]
            factors.append(convert_to_sympy_matrix(expr, expr.space))
            factors.extend([sympy.eye(s.dimension)
                            for s in all_spaces[own_space_index + 1:]])
            return tensor(*factors)
        if isinstance(expr, (Create, Jz, Jplus)):
            return SympyCreate(n)
        elif isinstance(expr, (Destroy, Jminus)):
            return SympyCreate(n).H
        elif isinstance(expr, Phase):
            phi = expr.phase
            result = sympy.zeros(n)
            for i in range(n):
                result[i, i] = sympy.exp(sympy.I * i * phi)
            return result
        elif isinstance(expr, Displace):
            alpha = expr.operands[1]
            a = SympyCreate(n)
            return (alpha * a - alpha.conjugate() * a.H).exp()
        elif isinstance(expr, Squeeze):
            eta = expr.operands[1]
            a = SympyCreate(n)
            return ((eta/2) * a**2 - (eta.conjugate()/2) * (a.H)**2).exp()
        elif isinstance(expr, LocalSigma):
            ket = basis_state(expr.index_j, n)
            bra = basis_state(expr.index_k, n).H
            return ket * bra
        else:
            raise ValueError("Cannot convert '%s' of type %s"
                             % (str(expr), type(expr)))
    elif (isinstance(expr, Operator) and isinstance(expr, Operation)):
        if isinstance(expr, OperatorPlus):
            s = convert_to_sympy_matrix(expr.operands[0], full_space)
            for op in expr.operands[1:]:
                s += convert_to_sympy_matrix(op, full_space)
            return s
        elif isinstance(expr, OperatorTimes):
            # if any factor acts non-locally, we need to expand distributively.
            if any(len(op.space) > 1 for op in expr.operands):
                se = expr.expand()
                if se == expr:
                    raise ValueError("Cannot represent as sympy matrix: %s"
                                     % expr)
                return convert_to_sympy_matrix(se, full_space)
            all_spaces = full_space.local_factors
            by_space = []
            ck = 0
            for ls in all_spaces:
                # group factors by associated local space
                ls_ops = [convert_to_sympy_matrix(o, o.space)
                          for o in expr.operands if o.space == ls]
                if len(ls_ops):
                    # compute factor associated with local space
                    by_space.append(ls_ops[0])
                    for ls_op in ls_ops[1:]:
                        by_space[-1] *= ls_op
                    ck += len(ls_ops)
                else:
                    # if trivial action, take identity matrix
                    by_space.append(sympy.eye(ls.dimension))
            assert ck == len(expr.operands)
            # combine local factors in tensor product
            if len(by_space) == 1:
                return by_space[0]
            else:
                return tensor(*by_space)

        elif isinstance(expr, Adjoint):