content: """
  This module contains functions for generating Qobj representation of a variety
  of commonly occuring quantum operators.
  """
  # Required for Sphinx to follow autodoc_type_aliases
  from __future__ import annotations
  
  __all__ = [
      'jmat', 'spin_Jx', 'spin_Jy', 'spin_Jz', 'spin_Jm', 'spin_Jp',
      'spin_J_set', 'sigmap', 'sigmam', 'sigmax', 'sigmay', 'sigmaz',
      'destroy', 'create', 'fdestroy', 'fcreate', 'qeye', 'identity',
      'position', 'momentum', 'num', 'squeeze', 'squeezing', 'displace',
      'commutator', 'qutrit_ops', 'qdiags', 'phase', 'qzero', 'charge',
      'tunneling', 'qft', 'qzero_like', 'qeye_like', 'swap',
  ]
  
  import numpy as np
  from typing import Literal, overload
  from . import data as _data
  from .qobj import Qobj
  from .dimensions import Space
  from .. import settings
  from ..typing import DimensionLike, SpaceLike, LayerType
  
  def qdiags(
      diagonals: np.typing.ArrayLike | list[np.typing.ArrayLike],
      offsets: int | list[int] = None,
      dims: DimensionLike = None,
      shape: tuple[int, int] = None,
      *,
      dtype: LayerType = None,
  ) -> Qobj:
      """
      Constructs an operator from an array of diagonals.
  
      Parameters
      ----------
      diagonals : array_like or sequence of array_like
          Array of elements to place along the selected diagonals.
  
      offsets : int or sequence of ints, optional
          Sequence for diagonals to be set:
              - k=0 main diagonal
              - k>0 kth upper diagonal
              - k<0 kth lower diagonal
  
      dims : list, optional
          Dimensions for operator
  
      shape : list, tuple, optional
          Shape of operator.  If omitted, a square operator large enough
          to contain the diagonals is generated.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Examples
      --------
      >>> qdiags(sqrt(range(1, 4)), 1) # doctest: +SKIP
      Quantum object: dims = [[4], [4]], \
  shape = [4, 4], type = oper, isherm = False
      Qobj data =
      [[ 0.          1.          0.          0.        ]
       [ 0.          0.          1.41421356  0.        ]
       [ 0.          0.          0.          1.73205081]
       [ 0.          0.          0.          0.        ]]
  
      """
      dtype = _data._parse_default_dtype(dtype, "diagonal")
      offsets = [0] if offsets is None else offsets
      if not isinstance(offsets, list):
          offsets = [offsets]
      if len(offsets) == 1 and offsets[0] != 0:
          isherm = False
          isunitary = False
      elif offsets == [0]:
          isherm = np.all(np.abs(np.imag(diagonals)) <= settings.core["atol"])
          isunitary = np.all(np.abs(diagonals) - 1 <= settings.core["atol"])
      else:
          isherm = None
          isunitary = None
      data = _data.diag[dtype](diagonals, offsets, shape)
      return Qobj(
          data, copy=False, dtype=dtype,
          dims=dims, isherm=isherm, isunitary=isunitary
      )
  
  
  @overload
  def jmat(
      j: float,
      which: Literal[None],
      *,
      dtype: LayerType = None
  ) -> tuple[Qobj]: ...
  
  @overload
  def jmat(
      j: float,
      which: Literal["x", "y", "z", "+", "-"],
      *,
      dtype: LayerType = None
  ) -> Qobj: ...
  
  def jmat(
      j: float,
      which: Literal["x", "y", "z", "+", "-", None] = None,
      *,
      dtype: LayerType = None
  ) -> Qobj | tuple[Qobj]:
      """Higher-order spin operators:
  
      Parameters
      ----------
      j : float
          Spin of operator
  
      which : str, optional
          Which operator to return 'x','y','z','+','-'.
          If not given, then output is ['x','y','z']
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      jmat : Qobj or tuple of Qobj
          ``qobj`` for requested spin operator(s).
  
  
      Examples
      --------
      >>> jmat(1) # doctest: +SKIP
      [ Quantum object: dims = [[3], [3]], \
  shape = [3, 3], type = oper, isHerm = True
      Qobj data =
      [[ 0.          0.70710678  0.        ]
       [ 0.70710678  0.          0.70710678]
       [ 0.          0.70710678  0.        ]]
       Quantum object: dims = [[3], [3]], \
  shape = [3, 3], type = oper, isHerm = True
      Qobj data =
      [[ 0.+0.j          0.-0.70710678j  0.+0.j        ]
       [ 0.+0.70710678j  0.+0.j          0.-0.70710678j]
       [ 0.+0.j          0.+0.70710678j  0.+0.j        ]]
       Quantum object: dims = [[3], [3]], \
  shape = [3, 3], type = oper, isHerm = True
      Qobj data =
      [[ 1.  0.  0.]
       [ 0.  0.  0.]
       [ 0.  0. -1.]]]
      """
      dtype = _data._parse_default_dtype(dtype, "sparse")
      if int(2 * j) != 2 * j or j < 0:
          raise ValueError('j must be a non-negative integer or half-integer')
  
      if not which:
          return (
              jmat(j, 'x', dtype=dtype),
              jmat(j, 'y', dtype=dtype),
              jmat(j, 'z', dtype=dtype)
          )
  
      dims = [[int(2*j + 1)]]*2
      if which == '+':
          return Qobj(_jplus(j, dtype=dtype), dims=dims, dtype=dtype,
                      isherm=False, isunitary=False, copy=False)
      if which == '-':
          return Qobj(_jplus(j, dtype=dtype).adjoint(), dims=dims, dtype=dtype,
                      isherm=False, isunitary=False, copy=False)
      if which == 'x':
          A = _jplus(j, dtype=dtype)
          return Qobj(_data.add(A, A.adjoint()) * 0.5, dims=dims, dtype=dtype,
                      isherm=True, isunitary=False, copy=False)
      if which == 'y':
          A = _data.mul(_jplus(j, dtype=dtype), -0.5j)
          return Qobj(_data.add(A, A.adjoint()), dims=dims, dtype=dtype,
                      isherm=True, isunitary=False, copy=False)
      if which == 'z':
          return Qobj(_jz(j, dtype=dtype), dims=dims, dtype=dtype,
                      isherm=True, isunitary=False, copy=False)
      raise ValueError('Invalid spin operator: ' + which)
  
  
  def _jplus(j, *, dtype=None):
      """
      Internal functions for generating the data representing the J-plus
      operator.
      """
      m = np.arange(j, -j - 1, -1, dtype=complex)
      data = np.sqrt(j * (j + 1) - m * (m + 1))[1:]
      return _data.diag[dtype](data, 1)
  
  
  def _jz(j, *, dtype=None):
      """
      Internal functions for generating the data representing the J-z operator.
      """
      dtype = _data._parse_default_dtype(dtype, "sparse")
      N = int(2*j + 1)
      data = np.array([j-k for k in range(N)], dtype=complex)
      return _data.diag[dtype](data, 0)
  
  
  #
  # Spin j operators:
  #
  def spin_Jx(j: float, *, dtype: LayerType = None) -> Qobj:
      """Spin-j x operator
  
      Parameters
      ----------
      j : float
          Spin of operator
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      op : Qobj
          ``qobj`` representation of the operator.
  
      """
      return jmat(j, 'x', dtype=dtype)
  
  
  def spin_Jy(j: float, *, dtype: LayerType = None) -> Qobj:
      """Spin-j y operator
  
      Parameters
      ----------
      j : float
          Spin of operator
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      op : Qobj
          ``qobj`` representation of the operator.
  
      """
      return jmat(j, 'y', dtype=dtype)
  
  
  def spin_Jz(j: float, *, dtype: LayerType = None) -> Qobj:
      """Spin-j z operator
  
      Parameters
      ----------
      j : float
          Spin of operator
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      op : Qobj
          ``qobj`` representation of the operator.
  
      """
      return jmat(j, 'z', dtype=dtype)
  
  
  def spin_Jm(j: float, *, dtype: LayerType = None) -> Qobj:
      """Spin-j annihilation operator
  
      Parameters
      ----------
      j : float
          Spin of operator
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      op : Qobj
          ``qobj`` representation of the operator.
  
      """
      return jmat(j, '-', dtype=dtype)
  
  
  def spin_Jp(j: float, *, dtype: LayerType = None) -> Qobj:
      """Spin-j creation operator
  
      Parameters
      ----------
      j : float
          Spin of operator
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      op : Qobj
          ``qobj`` representation of the operator.
  
      """
      return jmat(j, '+', dtype=dtype)
  
  
  def spin_J_set(j: float, *, dtype: LayerType = None) -> tuple[Qobj]:
      """Set of spin-j operators (x, y, z)
  
      Parameters
      ----------
      j : float
          Spin of operators
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      list : tuple of Qobj
          list of ``qobj`` representating of the spin operator.
  
      """
      return jmat(j, dtype=dtype)
  
  
  # Pauli spin-1/2 operators.
  #
  # These are so common in quantum information that we want them to be
  # near-instantaneous to initialise, so we cache them at package import, and
  # just return copies when someone requests one.
  _SIGMAP = jmat(0.5, '+')
  _SIGMAM = jmat(0.5, '-')
  _SIGMAX = 2 * jmat(0.5, 'x')
  _SIGMAX._isunitary = True
  _SIGMAY = 2 * jmat(0.5, 'y')
  _SIGMAY._isunitary = True
  _SIGMAZ = 2 * jmat(0.5, 'z')
  _SIGMAZ._isunitary = True
  
  
  def sigmap(*, dtype: LayerType = None) -> Qobj:
      """Creation operator for Pauli spins.
  
      Parameters
      ----------
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Examples
      --------
      >>> sigmap() # doctest: +SKIP
      Quantum object: dims = [[2], [2]], \
  shape = [2, 2], type = oper, isHerm = False
      Qobj data =
      [[ 0.  1.]
       [ 0.  0.]]
  
      """
      dtype = _data._parse_default_dtype(dtype, "sparse")
      return _SIGMAP.to(dtype, True)
  
  
  def sigmam(*, dtype: LayerType = None) -> Qobj:
      """Annihilation operator for Pauli spins.
  
      Parameters
      ----------
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Examples
      --------
      >>> sigmam() # doctest: +SKIP
      Quantum object: dims = [[2], [2]], \
  shape = [2, 2], type = oper, isHerm = False
      Qobj data =
      [[ 0.  0.]
       [ 1.  0.]]
  
      """
      dtype = _data._parse_default_dtype(dtype, "sparse")
      return _SIGMAM.to(dtype, True)
  
  
  def sigmax(*, dtype: LayerType = None) -> Qobj:
      """Pauli spin 1/2 sigma-x operator
  
      Parameters
      ----------
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Examples
      --------
      >>> sigmax() # doctest: +SKIP
      Quantum object: dims = [[2], [2]], \
  shape = [2, 2], type = oper, isHerm = False
      Qobj data =
      [[ 0.  1.]
       [ 1.  0.]]
  
      """
      dtype = _data._parse_default_dtype(dtype, "sparse")
      return _SIGMAX.to(dtype, True)
  
  
  def sigmay(*, dtype: LayerType = None) -> Qobj:
      """Pauli spin 1/2 sigma-y operator.
  
      Parameters
      ----------
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Examples
      --------
      >>> sigmay() # doctest: +SKIP
      Quantum object: dims = [[2], [2]], \
  shape = [2, 2], type = oper, isHerm = True
      Qobj data =
      [[ 0.+0.j  0.-1.j]
       [ 0.+1.j  0.+0.j]]
  
      """
      dtype = _data._parse_default_dtype(dtype, "sparse")
      return _SIGMAY.to(dtype, True)
  
  
  def sigmaz(*, dtype: LayerType = None) -> Qobj:
      """Pauli spin 1/2 sigma-z operator.
  
      Parameters
      ----------
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Examples
      --------
      >>> sigmaz() # doctest: +SKIP
      Quantum object: dims = [[2], [2]], \
  shape = [2, 2], type = oper, isHerm = True
      Qobj data =
      [[ 1.  0.]
       [ 0. -1.]]
  
      """
      dtype = _data._parse_default_dtype(dtype, "sparse")
      return _SIGMAZ.to(dtype, True)
  
  
  def destroy(N: int, offset: int = 0, *, dtype: LayerType = None) -> Qobj:
      """
      Destruction (lowering) operator.
  
      Parameters
      ----------
      N : int
          Number of basis states in the Hilbert space.
  
      offset : int, default: 0
          The lowest number state that is included in the finite number state
          representation of the operator.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper : qobj
          Qobj for lowering operator.
  
      Examples
      --------
      >>> destroy(4) # doctest: +SKIP
      Quantum object: dims=[[4], [4]], shape=(4, 4), type='oper', isherm=False
      Qobj data =
      [[ 0.00000000+0.j  1.00000000+0.j  0.00000000+0.j  0.00000000+0.j]
       [ 0.00000000+0.j  0.00000000+0.j  1.41421356+0.j  0.00000000+0.j]
       [ 0.00000000+0.j  0.00000000+0.j  0.00000000+0.j  1.73205081+0.j]
       [ 0.00000000+0.j  0.00000000+0.j  0.00000000+0.j  0.00000000+0.j]]
      """
      dtype = _data._parse_default_dtype(dtype, "diagonal")
      if not isinstance(N, (int, np.integer)):  # raise error if N not integer
          raise ValueError("Hilbert space dimension must be integer value")
      data = np.sqrt(np.arange(offset+1, N+offset, dtype=complex))
      return qdiags(data, 1, dtype=dtype)
  
  
  def create(N: int, offset: int = 0, *, dtype: LayerType = None) -> Qobj:
      """
      Creation (raising) operator.
  
      Parameters
      ----------
      N : int
          Number of basis states in the Hilbert space.
  
      offset : int, default: 0
          The lowest number state that is included in the finite number state
          representation of the operator.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper : qobj
          Qobj for raising operator.
  
      Examples
      --------
      >>> create(4) # doctest: +SKIP
      Quantum object: dims=[[4], [4]], shape=(4, 4), type='oper', isherm=False
      Qobj data =
      [[ 0.00000000+0.j  0.00000000+0.j  0.00000000+0.j  0.00000000+0.j]
       [ 1.00000000+0.j  0.00000000+0.j  0.00000000+0.j  0.00000000+0.j]
       [ 0.00000000+0.j  1.41421356+0.j  0.00000000+0.j  0.00000000+0.j]
       [ 0.00000000+0.j  0.00000000+0.j  1.73205081+0.j  0.00000000+0.j]]
      """
      dtype = _data._parse_default_dtype(dtype, "diagonal")
      if not isinstance(N, (int, np.integer)):  # raise error if N not integer
          raise ValueError("Hilbert space dimension must be integer value")
      data = np.sqrt(np.arange(offset+1, N+offset, dtype=complex))
      return qdiags(data, -1, dtype=dtype)
  
  
  def fdestroy(n_sites: int, site, dtype: LayerType = None) -> Qobj:
      """
      Fermionic destruction operator.
      We use the Jordan-Wigner transformation,
      making use of the Jordan-Wigner ZZ..Z strings,
      to construct this as follows:
  
      .. math::
  
          a_j = \\sigma_z^{\\otimes j} \\otimes
          (\\frac{\\sigma_x + i \\sigma_y}{2})
          \\otimes I^{\\otimes N-j-1}
  
      Parameters
      ----------
      n_sites : int
          Number of sites in Fock space.
  
      site : int
          The site in Fock space to add a fermion to.
          Corresponds to j in the above JW transform.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper : qobj
          Qobj for destruction operator.
  
      Examples
      --------
      >>> fdestroy(2) # doctest: +SKIP
      Quantum object: dims=[[2 2], [2 2]], shape=(4, 4), \
      type='oper', isherm=False
      Qobj data =
      [[0. 0. 1. 0.]
      [0. 0. 0. 1.]
      [0. 0. 0. 0.]
      [0. 0. 0. 0.]]
      """
      return _f_op(n_sites, site, 'destruction', dtype=dtype)
  
  
  def fcreate(n_sites: int, site, dtype: LayerType = None) -> Qobj:
      """
      Fermionic creation operator.
      We use the Jordan-Wigner transformation,
      making use of the Jordan-Wigner ZZ..Z strings,
      to construct this as follows:
  
      .. math::
  
          a_j = \\sigma_z^{\\otimes j}
          \\otimes (\\frac{\\sigma_x - i \\sigma_y}{2})
          \\otimes I^{\\otimes N-j-1}
  
  
      Parameters
      ----------
      n_sites : int
          Number of sites in Fock space.
  
      site : int
          The site in Fock space to add a fermion to.
          Corresponds to j in the above JW transform.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper : qobj
          Qobj for raising operator.
  
      Examples
      --------
      >>> fcreate(2) # doctest: +SKIP
      Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), \
      type = oper, isherm = False
      Qobj data =
      [[0. 0. 0. 0.]
      [0. 0. 0. 0.]
      [1. 0. 0. 0.]
      [0. 1. 0. 0.]]
      """
      return _f_op(n_sites, site, 'creation', dtype=dtype)
  
  
  def _f_op(n_sites, site, action, dtype: LayerType = None,):
      """ Makes fermionic creation and destruction operators.
      We use the Jordan-Wigner transformation,
      making use of the Jordan-Wigner ZZ..Z strings,
      to construct this as follows:
  
      .. math::
  
          a_j = \\sigma_z^{\\otimes j}
          \\otimes (frac{sigma_x \\pm i sigma_y}{2})
          \\otimes I^{\\otimes N-j-1}
  
      Parameters
      ----------
      action : str
          The type of operator to build.
          Can only be 'creation' or 'destruction'
  
      n_sites : int
          Number of sites in Fock space.
  
      site : int
          The site in Fock space to create/destroy a fermion on.
          Corresponds to j in the above JW transform.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper : qobj
          Qobj for destruction operator.
      """
      dtype = _data._parse_default_dtype(dtype, "sparse")
      # get `tensor` and sigma z objects
      from .tensor import tensor
      s_z = 2 * jmat(0.5, 'z', dtype=dtype)
  
      # sanity check
      if site < 0:
          raise ValueError(f'The specified site {site} cannot be \
                           less than 0.')
      elif 0 >= n_sites:
          raise ValueError(f'The specified number of sites {n_sites} \
                           cannot be equal to or less than 0.')
      elif site >= n_sites:
          raise ValueError(f'The specified site {site} is not in \
                           the range of {n_sites} sites.')
  
      # figure out which operator to build
      if action.lower() == 'creation':
          operator = create(2, dtype=dtype)
      elif action.lower() == 'destruction':
          operator = destroy(2, dtype=dtype)
      else:
          raise TypeError("Unknown operator '%s'. `action` must be \
                          either 'creation' or 'destruction.'" % action)
  
      eye = identity(2, dtype=dtype)
      opers = [s_z] * site + [operator] + [eye] * (n_sites - site - 1)
      out = tensor(opers).to(dtype)
      out.isherm = False
      out._isunitary = False
      return out
  
  
  def qzero(
      dimensions: SpaceLike,
      dims_right: SpaceLike = None,
      *,
      dtype: LayerType = None
  ) -> Qobj:
      """
      Zero operator.
  
      Parameters
      ----------
      dimensions : int, list of int, list of list of int, Space
          Number of basis states in the Hilbert space. If provided as a list of
          ints, then the dimension is the product over this list, but the
          ``dims`` property of the new Qobj are set to this list.  This can
          produce either `oper` or `super` depending on the passed `dimensions`.
  
      dims_right : int, list of int, list of list of int, Space, optional
          Number of basis states in the right Hilbert space when the operator is
          rectangular.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      qzero : qobj
          Zero operator Qobj.
  
      """
      dtype = _data._parse_default_dtype(dtype, "sparse")
      dims_left = Space(dimensions)
      size_left = dims_left.size
      if dims_right is None:
          dims_right = dims_left
          size_right = size_left
      else:
          dims_right = Space(dims_right)
          size_right = dims_right.size
      dims = [dims_left, dims_right]
      # A sparse matrix with no data is equal to a zero matrix.
      return Qobj(_data.zeros[dtype](size_left, size_right), dims=dims,
                  isherm=True, isunitary=False, copy=False, dtype=dtype)
  
  
  def qzero_like(qobj: Qobj) -> Qobj:
      """
      Zero operator of the same dims and type as the reference.
  
      Parameters
      ----------
      qobj : Qobj, QobjEvo
          Reference quantum object to copy the dims from.
  
      Returns
      -------
      qzero : qobj
          Zero operator Qobj.
  
      """
  
      return Qobj(
          _data.zeros[qobj.dtype](*qobj.shape), dims=qobj._dims,
          isherm=True, isunitary=False, copy=False, dtype=qobj.dtype
      )
  
  
  def qeye(dimensions: SpaceLike, *, dtype: LayerType = None) -> Qobj:
      """
      Identity operator.
  
      Parameters
      ----------
      dimensions : (int) or (list of int) or (list of list of int), Space
          Number of basis states in the Hilbert space. If provided as a list of
          ints, then the dimension is the product over this list, but the
          ``dims`` property of the new Qobj are set to this list.  This can
          produce either `oper` or `super` depending on the passed `dimensions`.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper : qobj
          Identity operator Qobj.
  
      Examples
      --------
      >>> qeye(3) # doctest: +SKIP
      Quantum object: dims = [[3], [3]], shape = (3, 3), type = oper, \
  isherm = True
      Qobj data =
      [[ 1.  0.  0.]
       [ 0.  1.  0.]
       [ 0.  0.  1.]]
      >>> qeye([2,2]) # doctest: +SKIP
      Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, \
  isherm = True
      Qobj data =
      [[1. 0. 0. 0.]
       [0. 1. 0. 0.]
       [0. 0. 1. 0.]
       [0. 0. 0. 1.]]
  
      """
      dtype = _data._parse_default_dtype(dtype, "diagonal")
      dimensions = Space(dimensions)
      return Qobj(_data.identity[dtype](dimensions.size), dims=[dimensions]*2,
                  isherm=True, isunitary=True, copy=False, dtype=dtype)
  
  
  # Name alias.
  identity = qeye
  
  
  def qeye_like(qobj: Qobj) -> Qobj:
      """
      Identity operator with the same dims and type as the reference quantum
      object.
  
      Parameters
      ----------
      qobj : Qobj, QobjEvo
          Reference quantum object to copy the dims from.
  
      Returns
      -------
      oper : qobj
          Identity operator Qobj.
  
      """
      if qobj.shape[0] != qobj.shape[1]:
          raise ValueError(
              "Can't create an identity matrix like a non square matrix."
          )
      return Qobj(
          _data.identity[qobj.dtype](qobj.shape[0]), dims=qobj._dims,
          isherm=True, isunitary=True, copy=False, dtype=qobj.dtype
      )
  
  
  def position(N: int, offset: int = 0, *, dtype: LayerType = None) -> Qobj:
      """
      Position operator :math:`x = 1 / sqrt(2) * (a + a.dag())`
  
      Parameters
      ----------
      N : int
          Number of basis states in the Hilbert space.
  
      offset : int, default: 0
          The lowest number state that is included in the finite number state
          representation of the operator.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper : qobj
          Position operator as Qobj.
      """
      dtype = _data._parse_default_dtype(dtype, "diagonal")
      a = destroy(N, offset=offset, dtype=dtype)
      position = np.sqrt(0.5) * (a + a.dag())
      position.isherm = True
      position._isunitary = False
      return position.to(dtype)
  
  
  def momentum(N: int, offset: int = 0, *, dtype: LayerType = None) -> Qobj:
      """
      Momentum operator p=-1j/sqrt(2)*(a-a.dag())
  
      Parameters
      ----------
      N : int
          Number of basis states in the Hilbert space.
  
      offset : int, default: 0
          The lowest number state that is included in the finite number state
          representation of the operator.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper : qobj
          Momentum operator as Qobj.
      """
      dtype = _data._parse_default_dtype(dtype, "diagonal")
      a = destroy(N, offset=offset, dtype=dtype)
      momentum = -1j * np.sqrt(0.5) * (a - a.dag())
      momentum.isherm = True
      momentum._isunitary = False
      return momentum.to(dtype)
  
  
  def num(N: int, offset: int = 0, *, dtype: LayerType = None) -> Qobj:
      """
      Quantum object for number operator.
  
      Parameters
      ----------
      N : int
          Number of basis states in the Hilbert space.
  
      offset : int, default: 0
          The lowest number state that is included in the finite number state
          representation of the operator.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper: qobj
          Qobj for number operator.
  
      Examples
      --------
      >>> num(4) # doctest: +SKIP
      Quantum object: dims=[[4], [4]], shape=(4, 4), type='oper', isherm=True
      Qobj data =
      [[0 0 0 0]
       [0 1 0 0]
       [0 0 2 0]
       [0 0 0 3]]
      """
      dtype = _data._parse_default_dtype(dtype, "diagonal")
      data = np.arange(offset, offset + N, dtype=complex)
      return qdiags(data, 0, dtype=dtype)
  
  
  def squeeze(
      N: int,
      z: float,
      offset: int = 0,
      *,
      dtype: LayerType = None,
  ) -> Qobj:
      """Single-mode squeezing operator.
  
      Parameters
      ----------
      N : int
          Dimension of hilbert space.
  
      z : float/complex
          Squeezing parameter.
  
      offset : int, default: 0
          The lowest number state that is included in the finite number state
          representation of the operator.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper : :class:`.Qobj`
          Squeezing operator.
  
  
      Examples
      --------
      >>> squeeze(4, 0.25) # doctest: +SKIP
      Quantum object: dims = [[4], [4]], \
  shape = [4, 4], type = oper, isHerm = False
      Qobj data =
      [[ 0.98441565+0.j  0.00000000+0.j  0.17585742+0.j  0.00000000+0.j]
       [ 0.00000000+0.j  0.95349007+0.j  0.00000000+0.j  0.30142443+0.j]
       [-0.17585742+0.j  0.00000000+0.j  0.98441565+0.j  0.00000000+0.j]
       [ 0.00000000+0.j -0.30142443+0.j  0.00000000+0.j  0.95349007+0.j]]
  
      """
      dtype = _data._parse_default_dtype(dtype, "dense")
      asq = destroy(N, offset=offset, dtype=dtype) ** 2
      op = 0.5*np.conj(z)*asq - 0.5*z*asq.dag()
      out = op.expm(dtype=dtype)
      out.isherm = (N == 2) or (z == 0.)
      out._isunitary = True
      return out
  
  
  def squeezing(a1: Qobj, a2: Qobj, z: float) -> Qobj:
      """Generalized squeezing operator.
  
      .. math::
  
          S(z) = \\exp\\left(\\frac{1}{2}\\left(z^*a_1a_2
          - za_1^\\dagger a_2^\\dagger\\right)\\right)
  
      Parameters
      ----------
      a1 : :class:`.Qobj`
          Operator 1.
  
      a2 : :class:`.Qobj`
          Operator 2.
  
      z : float/complex
          Squeezing parameter.
  
      Returns
      -------
      oper : :class:`.Qobj`
          Squeezing operator.
  
      """
      b = 0.5 * (np.conj(z)*(a1 @ a2) - z*(a1.dag() @ a2.dag()))
      return b.expm()
  
  
  def displace(
      N: int,
      alpha: float,
      offset: int = 0,
      *,
      dtype: LayerType = None,
  ) -> Qobj:
      """Single-mode displacement operator.
  
      Parameters
      ----------
      N : int
          Number of basis states in the Hilbert space.
  
      alpha : float/complex
          Displacement amplitude.
  
      offset : int, default: 0
          The lowest number state that is included in the finite number state
          representation of the operator.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper : qobj
          Displacement operator.
  
      Examples
      --------
      >>> displace(4,0.25) # doctest: +SKIP
      Quantum object: dims = [[4], [4]], \
  shape = [4, 4], type = oper, isHerm = False
      Qobj data =
      [[ 0.96923323+0.j -0.24230859+0.j  0.04282883+0.j -0.00626025+0.j]
       [ 0.24230859+0.j  0.90866411+0.j -0.33183303+0.j  0.07418172+0.j]
       [ 0.04282883+0.j  0.33183303+0.j  0.84809499+0.j -0.41083747+0.j]
       [ 0.00626025+0.j  0.07418172+0.j  0.41083747+0.j  0.90866411+0.j]]
  
      """
      dtype = _data._parse_default_dtype(dtype, "dense")
      a = destroy(N, offset=offset)
      out = (alpha * a.dag() - np.conj(alpha) * a).expm(dtype=dtype)
      out.isherm = (alpha == 0.)
      out._isunitary = True
      return out
  
  
  def commutator(
      A: Qobj,
      B: Qobj,
      kind: Literal["normal", "anti"] = "normal"
  ) -> Qobj:
      """
      Return the commutator of kind `kind` (normal, anti) of the
      two operators A and B.
  
      Parameters
      ----------
      A, B : :obj:`Qobj`, :obj:`QobjEvo`
          The operators to compute the commutator of.
  
      kind: str {"normal", "anti"}, default: "anti"
          Which kind of commutator to compute.
      """
      if kind == 'normal':
          return A @ B - B @ A
  
      elif kind == 'anti':
          return A @ B + B @ A
  
      else:
          raise TypeError("Unknown commutator kind '%s'" % kind)
  
  
  def qutrit_ops(*, dtype: LayerType = None) -> list[Qobj]:
      """
      Operators for a three level system (qutrit).
  
      Parameters
      ----------
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      opers: array
          `array` of qutrit operators.
  
      """
      from .states import qutrit_basis
  
      dtype = _data._parse_default_dtype(dtype, "sparse")
      out = []
      basis = qutrit_basis()
      for i in range(3):
          op = (basis[i] @ basis[i].dag()).to(dtype)
          op.isherm = True
          op._isunitary = False
          out.append(op)
      for i in range(3):
          op = (basis[i] @ basis[(i+1)%3].dag()).to(dtype)
          op.isherm = False
          op._isunitary = False
          out.append(op)
      return out
  
  
  def phase(N: int, phi0: float = 0, *, dtype: LayerType = None) -> Qobj:
      """
      Single-mode Pegg-Barnett phase operator.
  
      Parameters
      ----------
      N : int
          Number of basis states in the Hilbert space.
  
      phi0 : float, default: 0
          Reference phase.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      oper : qobj
          Phase operator with respect to reference phase.
  
      Notes
      -----
      The Pegg-Barnett phase operator is Hermitian on a truncated Hilbert space.
  
      """
      dtype = _data._parse_default_dtype(dtype, "dense")
      phim = phi0 + (2 * np.pi * np.arange(N)) / N  # discrete phase angles
      n = np.arange(N)[:, np.newaxis]
      states = np.array([np.sqrt(kk) / np.sqrt(N) * np.exp(1j * n * kk)
                         for kk in phim])
      ops = np.sum([np.outer(st, st.conj()) for st in states], axis=0)
      return Qobj(
          ops,
          isherm=True,
          isunitary=False,
          dims=[[N], [N]],
          copy=False
      ).to(dtype)
  
  
  def charge(
      Nmax: int,
      Nmin: int = None,
      frac: float = 1,
      *,
      dtype: LayerType = None
  ) -> Qobj:
      """
      Generate the diagonal charge operator over charge states
      from Nmin to Nmax.
  
      Parameters
      ----------
      Nmax : int
          Maximum charge state to consider.
  
      Nmin : int, default: -Nmax
          Lowest charge state to consider.
  
      frac : float, default: 1
          Specify fractional charge if needed.
  
      dtype : type or str, optional
          Storage representation. Any data-layer known to ``qutip.data.to`` is
          accepted.
  
      Returns
      -------
      C : Qobj
          Charge operator over [Nmin, Nmax].
  
      Notes
      -----
      .. versionadded:: 3.2
  
      """
      dtype = _data._parse_default_dtype(dtype, "diagonal")
      if Nmin is None:
          Nmin = -Nmax
      diag = frac * np.arange(Nmin, Nmax+1, dtype=float)
      out = qdiags(diag, 0, dtype=dtype)
      out._isunitary = (len(diag) <= 2) and np.all(np.abs(diag) == 1.)
      return out
  
  
  def tunneling(N: int, m: int = 1, *, dtype: LayerType = None) -> Qobj:
      r"""
      Tunneling operator with elements of the form
      :math:`\\sum |N><N+m| + |N+m><N|`.
  
      Parameters
      ----------
      N : int
  