# Install qiskit-terra=0.21 to use group_commuting()
# ! pip install qiskit-terra=0.21

from collections import defaultdict
from numbers import Number
from typing import Dict

from collections import defaultdict

import numpy as np
import retworkx as rx

from qiskit.exceptions import QiskitError
from qiskit.quantum_info.operators.custom_iterator import CustomIterator
from qiskit.quantum_info.operators.mixins import GroupMixin, LinearMixin
from qiskit.quantum_info.operators.symplectic.base_pauli import BasePauli
from qiskit.quantum_info.operators.symplectic.pauli import Pauli
from qiskit.quantum_info.operators.symplectic.pauli_table import PauliTable
from qiskit.quantum_info.operators.symplectic.stabilizer_table import StabilizerTable
from qiskit.quantum_info.operators.linear_op import LinearOp
from qiskit.utils.deprecation import deprecate_function


 



# if version incompativility occurs in terra 0.21, set PauliList and SparsePauliOp class below. 
class PauliList(BasePauli, LinearMixin, GroupMixin):
    r"""List of N-qubit Pauli operators.

    This class is an efficient representation of a list of
    :class:`Pauli` operators. It supports 1D numpy array indexing
    returning a :class:`Pauli` for integer indexes or a
    :class:`PauliList` for slice or list indices.

    **Initialization**

    A PauliList object can be initialized in several ways.

        ``PauliList(list[str])``
            where strings are same representation with :class:`~qiskit.quantum_info.Pauli`.

        ``PauliList(Pauli) and PauliList(list[Pauli])``
            where Pauli is :class:`~qiskit.quantum_info.Pauli`.

        ``PauliList.from_symplectic(z, x, phase)``
            where ``z`` and ``x`` are 2 dimensional boolean ``numpy.ndarrays`` and ``phase`` is
            an integer in ``[0, 1, 2, 3]``.

    For example,

    .. jupyter-execute::

        import numpy as np

        from qiskit.quantum_info import Pauli, PauliList

        # 1. init from list[str]
        pauli_list = PauliList(["II", "+ZI", "-iYY"])
        print("1. ", pauli_list)

        pauli1 = Pauli("iXI")
        pauli2 = Pauli("iZZ")

        # 2. init from Pauli
        print("2. ", PauliList(pauli1))

        # 3. init from list[Pauli]
        print("3. ", PauliList([pauli1, pauli2]))

        # 4. init from np.ndarray
        z = np.array([[True, True], [False, False]])
        x = np.array([[False, True], [True, False]])
        phase = np.array([0, 1])
        pauli_list = PauliList.from_symplectic(z, x, phase)
        print("4. ", pauli_list)

    **Data Access**

    The individual Paulis can be accessed and updated using the ``[]``
    operator which accepts integer, lists, or slices for selecting subsets
    of PauliList. If integer is given, it returns Pauli not PauliList.

    .. jupyter-execute::

        pauli_list = PauliList(["XX", "ZZ", "IZ"])
        print("Integer: ", repr(pauli_list[1]))
        print("List: ", repr(pauli_list[[0, 2]]))
        print("Slice: ", repr(pauli_list[0:2]))

    **Iteration**

    Rows in the Pauli table can be iterated over like a list. Iteration can
    also be done using the label or matrix representation of each row using the
    :meth:`label_iter` and :meth:`matrix_iter` methods.
    """

    # Set the max number of qubits * paulis before string truncation
    __truncate__ = 2000

    def __init__(self, data):
        """Initialize the PauliList.

        Args:
            data (Pauli or list): input data for Paulis. If input is a list each item in the list
                                  must be a Pauli object or Pauli str.

        Raises:
            QiskitError: if input array is invalid shape.

        Additional Information:
            The input array is not copied so multiple Pauli tables
            can share the same underlying array.
        """
        if isinstance(data, BasePauli):
            base_z, base_x, base_phase = data._z, data._x, data._phase
        elif isinstance(data, StabilizerTable):
            # Conversion from legacy StabilizerTable
            base_z, base_x, base_phase = self._from_array(data.Z, data.X, 2 * data.phase)
        elif isinstance(data, PauliTable):
            # Conversion from legacy PauliTable
            base_z, base_x, base_phase = self._from_array(data.Z, data.X)
        else:
            # Conversion as iterable of Paulis
            base_z, base_x, base_phase = self._from_paulis(data)

        # Initialize BasePauli
        super().__init__(base_z, base_x, base_phase)

    # ---------------------------------------------------------------------
    # Representation conversions
    # ---------------------------------------------------------------------

    @property
    def settings(self):
        """Return settings."""
        return {"data": self.to_labels()}

    def __array__(self, dtype=None):
        """Convert to numpy array"""
        # pylint: disable=unused-argument
        shape = (len(self),) + 2 * (2**self.num_qubits,)
        ret = np.zeros(shape, dtype=complex)
        for i, mat in enumerate(self.matrix_iter()):
            ret[i] = mat
        return ret

    @staticmethod
    def _from_paulis(data):
        """Construct a PauliList from a list of Pauli data.

        Args:
            data (iterable): list of Pauli data.

        Returns:
            PauliList: the constructed PauliList.

        Raises:
            QiskitError: If the input list is empty or contains invalid
            Pauli strings.
        """
        if not isinstance(data, (list, tuple, set, np.ndarray)):
            data = [data]
        num_paulis = len(data)
        if num_paulis == 0:
            raise QiskitError("Input Pauli list is empty.")
        paulis = []
        for i in data:
            if not isinstance(i, Pauli):
                paulis.append(Pauli(i))
            else:
                paulis.append(i)
        num_qubits = paulis[0].num_qubits
        base_z = np.zeros((num_paulis, num_qubits), dtype=bool)
        base_x = np.zeros((num_paulis, num_qubits), dtype=bool)
        base_phase = np.zeros(num_paulis, dtype=int)
        for i, pauli in enumerate(paulis):
            base_z[i] = pauli._z
            base_x[i] = pauli._x
            base_phase[i] = pauli._phase
        return base_z, base_x, base_phase

    def __repr__(self):
        """Display representation."""
        return self._truncated_str(True)

    def __str__(self):
        """Print representation."""
        return self._truncated_str(False)

    def _truncated_str(self, show_class):
        stop = self._num_paulis
        if self.__truncate__:
            max_paulis = self.__truncate__ // self.num_qubits
            if self._num_paulis > max_paulis:
                stop = max_paulis
        labels = [str(self[i]) for i in range(stop)]
        prefix = "PauliList(" if show_class else ""
        tail = ")" if show_class else ""
        if stop != self._num_paulis:
            suffix = ", ...]" + tail
        else:
            suffix = "]" + tail
        list_str = np.array2string(
            np.array(labels), threshold=stop + 1, separator=", ", prefix=prefix, suffix=suffix
        )
        return prefix + list_str[:-1] + suffix

    def __eq__(self, other):
        """Entrywise comparison of Pauli equality."""
        if not isinstance(other, PauliList):
            other = PauliList(other)
        if not isinstance(other, BasePauli):
            return False
        return self._eq(other)

    def equiv(self, other):
        """Entrywise comparison of Pauli equivalence up to global phase.

        Args:
            other (PauliList or Pauli): a comparison object.

        Returns:
            np.ndarray: An array of True or False for entrywise equivalence
                        of the current table.
        """
        if not isinstance(other, PauliList):
            other = PauliList(other)
        return np.all(self.z == other.z, axis=1) & np.all(self.x == other.x, axis=1)

    # ---------------------------------------------------------------------
    # Direct array access
    # ---------------------------------------------------------------------
    @property
    def phase(self):
        """Return the phase exponent of the PauliList."""
        # Convert internal ZX-phase convention to group phase convention
        return np.mod(self._phase - self._count_y(), 4)

    @phase.setter
    def phase(self, value):
        # Convert group phase convetion to internal ZX-phase convention
        self._phase[:] = np.mod(value + self._count_y(), 4)

    @property
    def x(self):
        """The x array for the symplectic representation."""
        return self._x

    @x.setter
    def x(self, val):
        self._x[:] = val

    @property
    def z(self):
        """The z array for the symplectic representation."""
        return self._z

    @z.setter
    def z(self, val):
        self._z[:] = val

    # ---------------------------------------------------------------------
    # Size Properties
    # ---------------------------------------------------------------------

    @property
    def shape(self):
        """The full shape of the :meth:`array`"""
        return self._num_paulis, self.num_qubits

    @property
    def size(self):
        """The number of Pauli rows in the table."""
        return self._num_paulis

    def __len__(self):
        """Return the number of Pauli rows in the table."""
        return self._num_paulis

    # ---------------------------------------------------------------------
    # Pauli Array methods
    # ---------------------------------------------------------------------

    def __getitem__(self, index):
        """Return a view of the PauliList."""
        # Returns a view of specified rows of the PauliList
        # This supports all slicing operations the underlying array supports.
        if isinstance(index, tuple):
            if len(index) == 1:
                index = index[0]
            elif len(index) > 2:
                raise IndexError(f"Invalid PauliList index {index}")

        # Row-only indexing
        if isinstance(index, (int, np.integer)):
            # Single Pauli
            return Pauli(
                BasePauli(
                    self._z[np.newaxis, index],
                    self._x[np.newaxis, index],
                    self._phase[np.newaxis, index],
                )
            )
        elif isinstance(index, (slice, list, np.ndarray)):
            # Sub-Table view
            return PauliList(BasePauli(self._z[index], self._x[index], self._phase[index]))

        # Row and Qubit indexing
        return PauliList((self._z[index], self._x[index], 0))

    def __setitem__(self, index, value):
        """Update PauliList."""
        if isinstance(index, tuple):
            if len(index) == 1:
                index = index[0]
            elif len(index) > 2:
                raise IndexError(f"Invalid PauliList index {index}")

        # Modify specified rows of the PauliList
        if not isinstance(value, PauliList):
            value = PauliList(value)

        self._z[index] = value._z
        self._x[index] = value._x
        if not isinstance(index, tuple):
            # Row-only indexing
            self._phase[index] = value._phase
        else:
            # Row and Qubit indexing
            self._phase[index[0]] += value._phase
            self._phase %= 4

    def delete(self, ind, qubit=False):
        """Return a copy with Pauli rows deleted from table.

        When deleting qubits the qubit index is the same as the
        column index of the underlying :attr:`X` and :attr:`Z` arrays.

        Args:
            ind (int or list): index(es) to delete.
            qubit (bool): if True delete qubit columns, otherwise delete
                          Pauli rows (Default: False).

        Returns:
            PauliList: the resulting table with the entries removed.

        Raises:
            QiskitError: if ind is out of bounds for the array size or
                         number of qubits.
        """
        if isinstance(ind, int):
            ind = [ind]

        # Row deletion
        if not qubit:
            if max(ind) >= len(self):
                raise QiskitError(
                    "Indices {} are not all less than the size"
                    " of the PauliList ({})".format(ind, len(self))
                )
            z = np.delete(self._z, ind, axis=0)
            x = np.delete(self._x, ind, axis=0)
            phase = np.delete(self._phase, ind)

            return PauliList(BasePauli(z, x, phase))

        # Column (qubit) deletion
        if max(ind) >= self.num_qubits:
            raise QiskitError(
                "Indices {} are not all less than the number of"
                " qubits in the PauliList ({})".format(ind, self.num_qubits)
            )
        z = np.delete(self._z, ind, axis=1)
        x = np.delete(self._x, ind, axis=1)
        # Use self.phase, not self._phase as deleting qubits can change the
        # ZX phase convention
        return PauliList.from_symplectic(z, x, self.phase)

    def insert(self, ind, value, qubit=False):
        """Insert Pauli's into the table.

        When inserting qubits the qubit index is the same as the
        column index of the underlying :attr:`X` and :attr:`Z` arrays.

        Args:
            ind (int): index to insert at.
            value (PauliList): values to insert.
            qubit (bool): if True delete qubit columns, otherwise delete
                          Pauli rows (Default: False).

        Returns:
            PauliList: the resulting table with the entries inserted.

        Raises:
            QiskitError: if the insertion index is invalid.
        """
        if not isinstance(ind, int):
            raise QiskitError("Insert index must be an integer.")

        if not isinstance(value, PauliList):
            value = PauliList(value)

        # Row insertion
        size = self._num_paulis
        if not qubit:
            if ind > size:
                raise QiskitError(
                    "Index {} is larger than the number of rows in the"
                    " PauliList ({}).".format(ind, size)
                )
            base_z = np.insert(self._z, ind, value._z, axis=0)
            base_x = np.insert(self._x, ind, value._x, axis=0)
            base_phase = np.insert(self._phase, ind, value._phase)
            return PauliList(BasePauli(base_z, base_x, base_phase))

        # Column insertion
        if ind > self.num_qubits:
            raise QiskitError(
                "Index {} is greater than number of qubits"
                " in the PauliList ({})".format(ind, self.num_qubits)
            )
        if len(value) == 1:
            # Pad blocks to correct size
            value_x = np.vstack(size * [value.x])
            value_z = np.vstack(size * [value.z])
            value_phase = np.vstack(size * [value.phase])
        elif len(value) == size:
            #  Blocks are already correct size
            value_x = value.x
            value_z = value.z
            value_phase = value.phase
        else:
            # Blocks are incorrect size
            raise QiskitError(
                "Input PauliList must have a single row, or"
                " the same number of rows as the Pauli Table"
                " ({}).".format(size)
            )
        # Build new array by blocks
        z = np.hstack([self.z[:, :ind], value_z, self.z[:, ind:]])
        x = np.hstack([self.x[:, :ind], value_x, self.x[:, ind:]])
        phase = self.phase + value_phase

        return PauliList.from_symplectic(z, x, phase)

    def argsort(self, weight=False, phase=False):
        """Return indices for sorting the rows of the table.

        The default sort method is lexicographic sorting by qubit number.
        By using the `weight` kwarg the output can additionally be sorted
        by the number of non-identity terms in the Pauli, where the set of
        all Pauli's of a given weight are still ordered lexicographically.

        Args:
            weight (bool): Optionally sort by weight if True (Default: False).
            phase (bool): Optionally sort by phase before weight or order
                          (Default: False).

        Returns:
            array: the indices for sorting the table.
        """
        # Get order of each Pauli using
        # I => 0, X => 1, Y => 2, Z => 3
        x = self.x
        z = self.z
        order = 1 * (x & ~z) + 2 * (x & z) + 3 * (~x & z)
        phases = self.phase
        # Optionally get the weight of Pauli
        # This is the number of non identity terms
        if weight:
            weights = np.sum(x | z, axis=1)

        # To preserve ordering between successive sorts we
        # are use the 'stable' sort method
        indices = np.arange(self._num_paulis)

        # Initial sort by phases
        sort_inds = phases.argsort(kind="stable")
        indices = indices[sort_inds]
        order = order[sort_inds]
        if phase:
            phases = phases[sort_inds]
        if weight:
            weights = weights[sort_inds]

        # Sort by order
        for i in range(self.num_qubits):
            sort_inds = order[:, i].argsort(kind="stable")
            order = order[sort_inds]
            indices = indices[sort_inds]
            if weight:
                weights = weights[sort_inds]
            if phase:
                phases = phases[sort_inds]

        # If using weights we implement a sort by total number
        # of non-identity Paulis
        if weight:
            sort_inds = weights.argsort(kind="stable")
            indices = indices[sort_inds]
            phases = phases[sort_inds]

        # If sorting by phase we perform a final sort by the phase value
        # of each pauli
        if phase:
            indices = indices[phases.argsort(kind="stable")]
        return indices

    def sort(self, weight=False, phase=False):
        """Sort the rows of the table.

        The default sort method is lexicographic sorting by qubit number.
        By using the `weight` kwarg the output can additionally be sorted
        by the number of non-identity terms in the Pauli, where the set of
        all Pauli's of a given weight are still ordered lexicographically.

        **Example**

        Consider sorting all a random ordering of all 2-qubit Paulis

        .. jupyter-execute::

            from numpy.random import shuffle
            from qiskit.quantum_info.operators import PauliList

            # 2-qubit labels
            labels = ['II', 'IX', 'IY', 'IZ', 'XI', 'XX', 'XY', 'XZ',
                      'YI', 'YX', 'YY', 'YZ', 'ZI', 'ZX', 'ZY', 'ZZ']
            # Shuffle Labels
            shuffle(labels)
            pt = PauliList(labels)
            print('Initial Ordering')
            print(pt)

            # Lexicographic Ordering
            srt = pt.sort()
            print('Lexicographically sorted')
            print(srt)

            # Weight Ordering
            srt = pt.sort(weight=True)
            print('Weight sorted')
            print(srt)

        Args:
            weight (bool): optionally sort by weight if True (Default: False).
            phase (bool): Optionally sort by phase before weight or order
                          (Default: False).

        Returns:
            PauliList: a sorted copy of the original table.
        """
        return self[self.argsort(weight=weight, phase=phase)]

    def unique(self, return_index=False, return_counts=False):
        """Return unique Paulis from the table.

        **Example**

        .. jupyter-execute::

            from qiskit.quantum_info.operators import PauliList

            pt = PauliList(['X', 'Y', '-X', 'I', 'I', 'Z', 'X', 'iZ'])
            unique = pt.unique()
            print(unique)

        Args:
            return_index (bool): If True, also return the indices that
                                 result in the unique array.
                                 (Default: False)
            return_counts (bool): If True, also return the number of times
                                  each unique item appears in the table.

        Returns:
            PauliList: unique
                the table of the unique rows.

            unique_indices: np.ndarray, optional
                The indices of the first occurrences of the unique values in
                the original array. Only provided if ``return_index`` is True.

            unique_counts: np.array, optional
                The number of times each of the unique values comes up in the
                original array. Only provided if ``return_counts`` is True.
        """
        # Check if we need to stack the phase array
        if np.any(self._phase != self._phase[0]):
            # Create a single array of Pauli's and phases for calling np.unique on
            # so that we treat different phased Pauli's as unique
            array = np.hstack([self._z, self._x, self.phase.reshape((self.phase.shape[0], 1))])
        else:
            # All Pauli's have the same phase so we only need to sort the array
            array = np.hstack([self._z, self._x])

        # Get indexes of unique entries
        if return_counts:
            _, index, counts = np.unique(array, return_index=True, return_counts=True, axis=0)
        else:
            _, index = np.unique(array, return_index=True, axis=0)

        # Sort the index so we return unique rows in the original array order
        sort_inds = index.argsort()
        index = index[sort_inds]
        unique = PauliList(BasePauli(self._z[index], self._x[index], self._phase[index]))

        # Concatinate return tuples
        ret = (unique,)
        if return_index:
            ret += (index,)
        if return_counts:
            ret += (counts[sort_inds],)
        if len(ret) == 1:
            return ret[0]
        return ret

    # ---------------------------------------------------------------------
    # BaseOperator methods
    # ---------------------------------------------------------------------

    def tensor(self, other):
        """Return the tensor product with each Pauli in the list.

        Args:
            other (PauliList): another PauliList.

        Returns:
            PauliList: the list of tensor product Paulis.

        Raises:
            QiskitError: if other cannot be converted to a PauliList, does
                         not have either 1 or the same number of Paulis as
                         the current list.
        """
        if not isinstance(other, PauliList):
            other = PauliList(other)
        return PauliList(super().tensor(other))

    def expand(self, other):
        """Return the expand product of each Pauli in the list.

        Args:
            other (PauliList): another PauliList.

        Returns:
            PauliList: the list of tensor product Paulis.

        Raises:
            QiskitError: if other cannot be converted to a PauliList, does
                         not have either 1 or the same number of Paulis as
                         the current list.
        """
        if not isinstance(other, PauliList):
            other = PauliList(other)
        if len(other) not in [1, len(self)]:
            raise QiskitError(
                "Incompatible PauliLists. Other list must "
                "have either 1 or the same number of Paulis."
            )
        return PauliList(super().expand(other))

    def compose(self, other, qargs=None, front=False, inplace=False):
        """Return the composition self∘other for each Pauli in the list.

        Args:
            other (PauliList): another PauliList.
            qargs (None or list): qubits to apply dot product on (Default: None).
            front (bool): If True use `dot` composition method [default: False].
            inplace (bool): If True update in-place (default: False).

        Returns:
            PauliList: the list of composed Paulis.

        Raises:
            QiskitError: if other cannot be converted to a PauliList, does
                         not have either 1 or the same number of Paulis as
                         the current list, or has the wrong number of qubits
                         for the specified qargs.
        """
        if qargs is None:
            qargs = getattr(other, "qargs", None)
        if not isinstance(other, PauliList):
            other = PauliList(other)
        if len(other) not in [1, len(self)]:
            raise QiskitError(
                "Incompatible PauliLists. Other list must "
                "have either 1 or the same number of Paulis."
            )
        return PauliList(super().compose(other, qargs=qargs, front=front, inplace=inplace))

    # pylint: disable=arguments-differ
    def dot(self, other, qargs=None, inplace=False):
        """Return the composition other∘self for each Pauli in the list.

        Args:
            other (PauliList): another PauliList.
            qargs (None or list): qubits to apply dot product on (Default: None).
            inplace (bool): If True update in-place (default: False).

        Returns:
            PauliList: the list of composed Paulis.

        Raises:
            QiskitError: if other cannot be converted to a PauliList, does
                         not have either 1 or the same number of Paulis as
                         the current list, or has the wrong number of qubits
                         for the specified qargs.
        """
        return self.compose(other, qargs=qargs, front=True, inplace=inplace)

    def _add(self, other, qargs=None):
        """Append two PauliLists.

        If ``qargs`` are specified the other operator will be added
        assuming it is identity on all other subsystems.

        Args:
            other (PauliList): another table.
            qargs (None or list): optional subsystems to add on
                                  (Default: None)

        Returns:
            PauliList: the concatenated list self + other.
        """
        if qargs is None:
            qargs = getattr(other, "qargs", None)

        if not isinstance(other, PauliList):
            other = PauliList(other)

        self._op_shape._validate_add(other._op_shape, qargs)

        base_phase = np.hstack((self._phase, other._phase))

        if qargs is None or (sorted(qargs) == qargs and len(qargs) == self.num_qubits):
            base_z = np.vstack([self._z, other._z])
            base_x = np.vstack([self._x, other._x])
        else:
            # Pad other with identity and then add
            padded = BasePauli(
                np.zeros((other.size, self.num_qubits), dtype=bool),
                np.zeros((other.size, self.num_qubits), dtype=bool),
                np.zeros(other.size, dtype=int),
            )
            padded = padded.compose(other, qargs=qargs, inplace=True)
            base_z = np.vstack([self._z, padded._z])
            base_x = np.vstack([self._x, padded._x])

        return PauliList(BasePauli(base_z, base_x, base_phase))

    def _multiply(self, other):
        """Multiply each Pauli in the list by a phase.

        Args:
            other (complex or array): a complex number in [1, -1j, -1, 1j]

        Returns:
            PauliList: the list of Paulis other * self.

        Raises:
            QiskitError: if the phase is not in the set [1, -1j, -1, 1j].
        """
        return PauliList(super()._multiply(other))

    def conjugate(self):
        """Return the conjugate of each Pauli in the list."""
        return PauliList(super().conjugate())

    def transpose(self):
        """Return the transpose of each Pauli in the list."""
        return PauliList(super().transpose())

    def adjoint(self):
        """Return the adjoint of each Pauli in the list."""
        return PauliList(super().adjoint())

    def inverse(self):
        """Return the inverse of each Pauli in the list."""
        return PauliList(super().adjoint())

    # ---------------------------------------------------------------------
    # Utility methods
    # ---------------------------------------------------------------------

    def commutes(self, other, qargs=None):
        """Return True for each Pauli that commutes with other.

        Args:
            other (PauliList): another PauliList operator.
            qargs (list): qubits to apply dot product on (default: None).

        Returns:
            bool: True if Pauli's commute, False if they anti-commute.
        """
        if qargs is None:
            qargs = getattr(other, "qargs", None)
        if not isinstance(other, BasePauli):
            other = PauliList(other)
        return super().commutes(other, qargs=qargs)

    def anticommutes(self, other, qargs=None):
        """Return True if other Pauli that anticommutes with other.

        Args:
            other (PauliList): another PauliList operator.
            qargs (list): qubits to apply dot product on (default: None).

        Returns:
            bool: True if Pauli's anticommute, False if they commute.
        """
        return np.logical_not(self.commutes(other, qargs=qargs))

    def commutes_with_all(self, other):
        """Return indexes of rows that commute other.

        If other is a multi-row Pauli list the returned vector indexes rows
        of the current PauliList that commute with *all* Pauli's in other.
        If no rows satisfy the condition the returned array will be empty.

        Args:
            other (PauliList): a single Pauli or multi-row PauliList.

        Returns:
            array: index array of the commuting rows.
        """
        return self._commutes_with_all(other)

    def anticommutes_with_all(self, other):
        """Return indexes of rows that commute other.

        If other is a multi-row Pauli list the returned vector indexes rows
        of the current PauliList that anti-commute with *all* Pauli's in other.
        If no rows satisfy the condition the returned array will be empty.

        Args:
            other (PauliList): a single Pauli or multi-row PauliList.

        Returns:
            array: index array of the anti-commuting rows.
        """
        return self._commutes_with_all(other, anti=True)

    def _commutes_with_all(self, other, anti=False):
        """Return row indexes that commute with all rows in another PauliList.

        Args:
            other (PauliList): a PauliList.
            anti (bool): if True return rows that anti-commute, otherwise
                         return rows that commute (Default: False).

        Returns:
            array: index array of commuting or anti-commuting row.
        """
        if not isinstance(other, PauliList):
            other = PauliList(other)
        comms = self.commutes(other[0])
        (inds,) = np.where(comms == int(not anti))
        for pauli in other[1:]:
            comms = self[inds].commutes(pauli)
            (new_inds,) = np.where(comms == int(not anti))
            if new_inds.size == 0:
                # No commuting rows
                return new_inds
            inds = inds[new_inds]
        return inds

    def evolve(self, other, qargs=None, frame="h"):
        r"""Evolve the Pauli by a Clifford.

        This returns the Pauli :math:`P^\prime = C.P.C^\dagger`.

        By choosing the parameter frame='s', this function returns the Schrödinger evolution of the Pauli
        :math:`P^\prime = C.P.C^\dagger`. This option yields a faster calculation.

        Args:
            other (Pauli or Clifford or QuantumCircuit): The Clifford operator to evolve by.
            qargs (list): a list of qubits to apply the Clifford to.
            frame (string): 'h' for Heisenberg or 's' for Schrödinger framework.

        Returns:
            Pauli: the Pauli :math:`C.P.C^\dagger`.

        Raises:
            QiskitError: if the Clifford number of qubits and qargs don't match.
        """
        from qiskit.circuit import Instruction, QuantumCircuit
        from qiskit.quantum_info.operators.symplectic.clifford import Clifford

        if qargs is None:
            qargs = getattr(other, "qargs", None)

        if not isinstance(other, (BasePauli, Instruction, QuantumCircuit, Clifford)):
            # Convert to a PauliList
            other = PauliList(other)

        return PauliList(super().evolve(other, qargs=qargs, frame=frame))

    def to_labels(self, array=False):
        r"""Convert a PauliList to a list Pauli string labels.

        For large PauliLists converting using the ``array=True``
        kwarg will be more efficient since it allocates memory for
        the full Numpy array of labels in advance.

        .. list-table:: Pauli Representations
            :header-rows: 1

            * - Label
              - Symplectic
              - Matrix
            * - ``"I"``
              - :math:`[0, 0]`
              - :math:`\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}`
            * - ``"X"``
              - :math:`[1, 0]`
              - :math:`\begin{bmatrix} 0 & 1 \\ 1 & 0  \end{bmatrix}`
            * - ``"Y"``
              - :math:`[1, 1]`
              - :math:`\begin{bmatrix} 0 & -i \\ i & 0  \end{bmatrix}`
            * - ``"Z"``
              - :math:`[0, 1]`
              - :math:`\begin{bmatrix} 1 & 0 \\ 0 & -1  \end{bmatrix}`

        Args:
            array (bool): return a Numpy array if True, otherwise
                          return a list (Default: False).

        Returns:
            list or array: The rows of the PauliList in label form.
        """
        if (self.phase == 1).any():
            prefix_len = 2
        elif (self.phase > 0).any():
            prefix_len = 1
        else:
            prefix_len = 0
        str_len = self.num_qubits + prefix_len
        ret = np.zeros(self.size, dtype=f"<U{str_len}")
        iterator = self.label_iter()
        for i in range(self.size):
            ret[i] = next(iterator)
        if array:
            return ret
        return ret.tolist()

    def to_matrix(self, sparse=False, array=False):
        r"""Convert to a list or array of Pauli matrices.

        For large PauliLists converting using the ``array=True``
        kwarg will be more efficient since it allocates memory a full
        rank-3 Numpy array of matrices in advance.

        .. list-table:: Pauli Representations
            :header-rows: 1

            * - Label
              - Symplectic
              - Matrix
            * - ``"I"``
              - :math:`[0, 0]`
              - :math:`\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}`
            * - ``"X"``
              - :math:`[1, 0]`
              - :math:`\begin{bmatrix} 0 & 1 \\ 1 & 0  \end{bmatrix}`
            * - ``"Y"``
              - :math:`[1, 1]`
              - :math:`\begin{bmatrix} 0 & -i \\ i & 0  \end{bmatrix}`
            * - ``"Z"``
              - :math:`[0, 1]`
              - :math:`\begin{bmatrix} 1 & 0 \\ 0 & -1  \end{bmatrix}`

        Args:
            sparse (bool): if True return sparse CSR matrices, otherwise
                           return dense Numpy arrays (Default: False).
            array (bool): return as rank-3 numpy array if True, otherwise
                          return a list of Numpy arrays (Default: False).

        Returns:
            list: A list of dense Pauli matrices if `array=False` and `sparse=False`.
            list: A list of sparse Pauli matrices if `array=False` and `sparse=True`.
            array: A dense rank-3 array of Pauli matrices if `array=True`.
        """
        if not array:
            # We return a list of Numpy array matrices
            return list(self.matrix_iter(sparse=sparse))
        # For efficiency we also allow returning a single rank-3
        # array where first index is the Pauli row, and second two
        # indices are the matrix indices
        dim = 2**self.num_qubits
        ret = np.zeros((self.size, dim, dim), dtype=complex)
        iterator = self.matrix_iter(sparse=sparse)
        for i in range(self.size):
            ret[i] = next(iterator)
        return ret

    # ---------------------------------------------------------------------
    # Custom Iterators
    # ---------------------------------------------------------------------

    def label_iter(self):
        """Return a label representation iterator.

        This is a lazy iterator that converts each row into the string
        label only as it is used. To convert the entire table to labels use
        the :meth:`to_labels` method.

        Returns:
            LabelIterator: label iterator object for the PauliList.
        """

        class LabelIterator(CustomIterator):
            """Label representation iteration and item access."""

            def __repr__(self):
                return f"<PauliList_label_iterator at {hex(id(self))}>"

            def __getitem__(self, key):
                return self.obj._to_label(self.obj._z[key], self.obj._x[key], self.obj._phase[key])

        return LabelIterator(self)

    def matrix_iter(self, sparse=False):
        """Return a matrix representation iterator.

        This is a lazy iterator that converts each row into the Pauli matrix
        representation only as it is used. To convert the entire table to
        matrices use the :meth:`to_matrix` method.

        Args:
            sparse (bool): optionally return sparse CSR matrices if True,
                           otherwise return Numpy array matrices
                           (Default: False)

        Returns:
            MatrixIterator: matrix iterator object for the PauliList.
        """

        class MatrixIterator(CustomIterator):
            """Matrix representation iteration and item access."""

            def __repr__(self):
                return f"<PauliList_matrix_iterator at {hex(id(self))}>"

            def __getitem__(self, key):
                return self.obj._to_matrix(
                    self.obj._z[key], self.obj._x[key], self.obj._phase[key], sparse=sparse
                )

        return MatrixIterator(self)

    # ---------------------------------------------------------------------
    # Class methods
    # ---------------------------------------------------------------------

    @classmethod
    def from_symplectic(cls, z, x, phase=0):
        """Construct a PauliList from a symplectic data.

        Args:
            z (np.ndarray): 2D boolean Numpy array.
            x (np.ndarray): 2D boolean Numpy array.
            phase (np.ndarray or None): Optional, 1D integer array from Z_4.

        Returns:
            PauliList: the constructed PauliList.
        """
        base_z, base_x, base_phase = cls._from_array(z, x, phase)
        return cls(BasePauli(base_z, base_x, base_phase))

    def _noncommutation_graph(self, qubit_wise):
        """Create an edge list representing the non-commutation graph (Pauli Graph).

        An edge (i, j) is present if i and j are not commutable.

        Args:
            qubit_wise (bool): whether the commutation rule is applied to the whole operator,
                or on a per-qubit basis.

        Returns:
            List[Tuple(int,int)]: A list of pairs of indices of the PauliList that are not commutable.
        """
        # convert a Pauli operator into int vector where {I: 0, X: 2, Y: 3, Z: 1}
        mat1 = np.array(
            [op.z + 2 * op.x for op in self],
            dtype=np.int8,
        )
        mat2 = mat1[:, None]
        qubit_commutation_mat = (mat1 * mat2) * (mat1 - mat2)
        # adjacency_mat[i, j] is True if an edge (i, j) is present, i.e. i and j are non-commuting
        if qubit_wise:
            adjacency_mat = np.logical_or.reduce(qubit_commutation_mat, axis=2)
        else:
            adjacency_mat = np.logical_xor.reduce(qubit_commutation_mat, axis=2)
        # convert into list where tuple elements are non-commuting operators
        return list(zip(*np.where(np.triu(adjacency_mat, k=1))))

    def _create_graph(self, qubit_wise):
        """Transform measurement operator grouping problem into graph coloring problem

        Args:
            qubit_wise (bool): whether the commutation rule is applied to the whole operator,
                or on a per-qubit basis.

        Returns:
            retworkx.PyGraph: A class of undirected graphs
        """

        edges = self._noncommutation_graph(qubit_wise)
        graph = rx.PyGraph()
        graph.add_nodes_from(range(self.size))
        graph.add_edges_from_no_data(edges)
        return graph

    def group_qubit_wise_commuting(self):
        """Partition a PauliList into sets of mutually qubit-wise commuting Pauli strings.

        Returns:
            List[PauliList]: List of PauliLists where each PauliList contains commutable Pauli operators.
        """
        return self.group_commuting(qubit_wise=True)

    def group_commuting(self, qubit_wise=False):
        """Partition a PauliList into sets of commuting Pauli strings.

        Args:
            qubit_wise (bool): whether the commutation rule is applied to the whole operator,
                or on a per-qubit basis.  For example:

                .. code-block:: python

                    >>> op = PauliList(["XX", "YY", "IZ", "ZZ"])
                    >>> op.group_commuting()
                    [PauliList(['XX', 'YY']), PauliList(['IZ', 'ZZ'])]
                    >>> op.group_commuting(qubit_wise=True)
                    [PauliList(['XX']), PauliList(['YY']), PauliList(['IZ', 'ZZ'])]

        Returns:
            List[PauliList]: List of PauliLists where each PauliList contains commuting Pauli operators.
        """

        graph = self._create_graph(qubit_wise)
        # Keys in coloring_dict are nodes, values are colors
        coloring_dict = rx.graph_greedy_color(graph)
        groups = defaultdict(list)
        for idx, color in coloring_dict.items():
            groups[color].append(idx)
        return [self[group] for group in groups.values()]
    
class SparsePauliOp(LinearOp):
    """Sparse N-qubit operator in a Pauli basis representation.

    This is a sparse representation of an N-qubit matrix
    :class:`~qiskit.quantum_info.Operator` in terms of N-qubit
    :class:`~qiskit.quantum_info.PauliList` and complex coefficients.

    It can be used for performing operator arithmetic for hundred of qubits
    if the number of non-zero Pauli basis terms is sufficiently small.

    The Pauli basis components are stored as a
    :class:`~qiskit.quantum_info.PauliList` object and can be accessed
    using the :attr:`~SparsePauliOp.paulis` attribute. The coefficients
    are stored as a complex Numpy array vector and can be accessed using
    the :attr:`~SparsePauliOp.coeffs` attribute.
    """

    def __init__(self, data, coeffs=None, *, ignore_pauli_phase=False, copy=True):
        """Initialize an operator object.

        Args:
            data (PauliList or SparsePauliOp or PauliTable or Pauli or list or str): Pauli list of
                terms.  A list of Pauli strings or a Pauli string is also allowed.
            coeffs (np.ndarray): complex coefficients for Pauli terms.

                .. note::

                    If ``data`` is a :obj:`~SparsePauliOp` and ``coeffs`` is not ``None``, the value
                    of the ``SparsePauliOp.coeffs`` will be ignored, and only the passed keyword
                    argument ``coeffs`` will be used.

            ignore_pauli_phase (bool): if true, any ``phase`` component of a given :obj:`~PauliList`
                will be assumed to be zero.  This is more efficient in cases where a
                :obj:`~PauliList` has been constructed purely for this object, and it is already
                known that the phases in the ZX-convention are zero.  It only makes sense to pass
                this option when giving :obj:`~PauliList` data.  (Default: False)
            copy (bool): copy the input data if True, otherwise assign it directly, if possible.
                (Default: True)

        Raises:
            QiskitError: If the input data or coeffs are invalid.
        """
        if ignore_pauli_phase and not isinstance(data, PauliList):
            raise QiskitError("ignore_pauli_list=True is only valid with PauliList data")

        if isinstance(data, SparsePauliOp):
            if coeffs is None:
                coeffs = data.coeffs
            data = data._pauli_list
            # `SparsePauliOp._pauli_list` is already compatible with the internal ZX-phase
            # convention.  See `BasePauli._from_array` for the internal ZX-phase convention.
            ignore_pauli_phase = True

        pauli_list = PauliList(data.copy() if copy and hasattr(data, "copy") else data)

        if coeffs is None:
            coeffs = np.ones(pauli_list.size, dtype=complex)
        else:
            coeffs = np.array(coeffs, copy=copy, dtype=complex)

        if ignore_pauli_phase:
            # Fast path used in copy operations, where the phase of the PauliList is already known
            # to be zero.  This path only works if the input data is compatible with the internal
            # ZX-phase convention.
            self._pauli_list = pauli_list
            self._coeffs = coeffs
        else:
            # move the phase of `pauli_list` to `self._coeffs`
            phase = pauli_list.phase
            self._coeffs = np.asarray((-1j) ** phase * coeffs, dtype=complex)
            pauli_list._phase = np.mod(pauli_list._phase - phase, 4)
            self._pauli_list = pauli_list

        if self._coeffs.shape != (self._pauli_list.size,):
            raise QiskitError(
                "coeff vector is incorrect shape for number"
                " of Paulis {} != {}".format(self._coeffs.shape, self._pauli_list.size)
            )
        # Initialize LinearOp
        super().__init__(num_qubits=self._pauli_list.num_qubits)

    def __array__(self, dtype=None):
        if dtype:
            return np.asarray(self.to_matrix(), dtype=dtype)
        return self.to_matrix()

    def __repr__(self):
        prefix = "SparsePauliOp("
        pad = len(prefix) * " "
        return "{}{},\n{}coeffs={})".format(
            prefix,
            self.paulis.to_labels(),
            pad,
            np.array2string(self.coeffs, separator=", "),
        )

    def __eq__(self, other):
        """Entrywise comparison of two SparsePauliOp operators"""
        return (
            super().__eq__(other)
            and self.coeffs.shape == other.coeffs.shape
            and np.allclose(self.coeffs, other.coeffs)
            and self.paulis == other.paulis
        )

    def equiv(self, other):
        """Check if two SparsePauliOp operators are equivalent.

        Args:
            other (SparsePauliOp): an operator object.

        Returns:
            bool: True if the operator is equivalent to ``self``.
        """
        if not super().__eq__(other):
            return False
        return np.allclose((self - other).simplify().coeffs, [0])

    @property
    def settings(self) -> Dict:
        """Return settings."""
        return {"data": self._pauli_list, "coeffs": self._coeffs}

    # ---------------------------------------------------------------------
    # Data accessors
    # ---------------------------------------------------------------------

    @property
    def size(self):
        """The number of Pauli of Pauli terms in the operator."""
        return self._pauli_list.size

    def __len__(self):
        """Return the size."""
        return self.size

    # pylint: disable=bad-docstring-quotes

    @property
    @deprecate_function(
        "The SparsePauliOp.table method is deprecated as of Qiskit Terra 0.19.0 "
        "and will be removed no sooner than 3 months after the releasedate. "
        "Use SparsePauliOp.paulis method instead.",
    )
    def table(self):
        """DEPRECATED - Return the the PauliTable."""
        return PauliTable(np.column_stack((self.paulis.x, self.paulis.z)))

    @table.setter
    @deprecate_function(
        "The SparsePauliOp.table method is deprecated as of Qiskit Terra 0.19.0 "
        "and will be removed no sooner than 3 months after the releasedate. "
        "Use SparsePauliOp.paulis method instead.",
    )
    def table(self, value):
        if not isinstance(value, PauliTable):
            value = PauliTable(value)
        self._pauli_list = PauliList(value)

    @property
    def paulis(self):
        """Return the the PauliList."""
        return self._pauli_list

    @paulis.setter
    def paulis(self, value):
        if not isinstance(value, PauliList):
            value = PauliList(value)
        self._pauli_list = value

    @property
    def coeffs(self):
        """Return the Pauli coefficients."""
        return self._coeffs

    @coeffs.setter
    def coeffs(self, value):
        """Set Pauli coefficients."""
        self._coeffs[:] = value

    def __getitem__(self, key):
        """Return a view of the SparsePauliOp."""
        # Returns a view of specified rows of the PauliList
        # This supports all slicing operations the underlying array supports.
        if isinstance(key, (int, np.integer)):
            key = [key]
        return SparsePauliOp(self.paulis[key], self.coeffs[key])

    def __setitem__(self, key, value):
        """Update SparsePauliOp."""
        # Modify specified rows of the PauliList
        if not isinstance(value, SparsePauliOp):
            value = SparsePauliOp(value)
        self.paulis[key] = value.paulis
        self.coeffs[key] = value.coeffs

    # ---------------------------------------------------------------------
    # LinearOp Methods
    # ---------------------------------------------------------------------

    def conjugate(self):
        # Conjugation conjugates phases and also Y.conj() = -Y
        # Hence we need to multiply conjugated coeffs by -1
        # for rows with an odd number of Y terms.
        # Find rows with Ys
        ret = self.transpose()
        ret._coeffs = ret._coeffs.conj()
        return ret

    def transpose(self):
        # The only effect transposition has is Y.T = -Y
        # Hence we need to multiply coeffs by -1 for rows with an
        # odd number of Y terms.
        ret = self.copy()
        minus = (-1) ** np.mod(np.sum(ret.paulis.x & ret.paulis.z, axis=1), 2)
        ret._coeffs *= minus
        return ret

    def adjoint(self):
        # Pauli's are self adjoint, so we only need to
        # conjugate the phases
        ret = self.copy()
        ret._coeffs = ret._coeffs.conj()
        return ret

    def compose(self, other, qargs=None, front=False):
        if qargs is None:
            qargs = getattr(other, "qargs", None)

        if not isinstance(other, SparsePauliOp):
            other = SparsePauliOp(other)

        # Validate composition dimensions and qargs match
        self._op_shape.compose(other._op_shape, qargs, front)

        if qargs is not None:
            x1, z1 = self.paulis.x[:, qargs], self.paulis.z[:, qargs]
        else:
            x1, z1 = self.paulis.x, self.paulis.z
        x2, z2 = other.paulis.x, other.paulis.z
        num_qubits = other.num_qubits

        # This method is the outer version of `BasePauli.compose`.
        # `x1` and `z1` have shape `(self.size, num_qubits)`.
        # `x2` and `z2` have shape `(other.size, num_qubits)`.
        # `x1[:, no.newaxis]` results in shape `(self.size, 1, num_qubits)`.
        # `ar = ufunc(x1[:, np.newaxis], x2)` will be in shape `(self.size, other.size, num_qubits)`.
        # So, `ar.reshape((-1, num_qubits))` will be in shape `(self.size * other.size, num_qubits)`.
        # Ref: https://numpy.org/doc/stable/user/theory.broadcasting.html

        phase = np.add.outer(self.paulis._phase, other.paulis._phase).reshape(-1)
        if front:
            q = np.logical_and(x1[:, np.newaxis], z2).reshape((-1, num_qubits))
        else:
            q = np.logical_and(z1[:, np.newaxis], x2).reshape((-1, num_qubits))
        # `np.mod` will be applied to `phase` in `SparsePauliOp.__init__`
        phase = phase + 2 * np.sum(q, axis=1)

        x3 = np.logical_xor(x1[:, np.newaxis], x2).reshape((-1, num_qubits))
        z3 = np.logical_xor(z1[:, np.newaxis], z2).reshape((-1, num_qubits))

        if qargs is None:
            pauli_list = PauliList(BasePauli(z3, x3, phase))
        else:
            x4 = np.repeat(self.paulis.x, other.size, axis=0)
            z4 = np.repeat(self.paulis.z, other.size, axis=0)
            x4[:, qargs] = x3
            z4[:, qargs] = z3
            pauli_list = PauliList(BasePauli(z4, x4, phase))

        # note: the following is a faster code equivalent to
        # `coeffs = np.kron(self.coeffs, other.coeffs)`
        # since `self.coeffs` and `other.coeffs` are both 1d arrays.
        coeffs = np.multiply.outer(self.coeffs, other.coeffs).ravel()
        return SparsePauliOp(pauli_list, coeffs, copy=False)

    def tensor(self, other):
        if not isinstance(other, SparsePauliOp):
            other = SparsePauliOp(other)
        return self._tensor(self, other)

    def expand(self, other):
        if not isinstance(other, SparsePauliOp):
            other = SparsePauliOp(other)
        return self._tensor(other, self)

    @classmethod
    def _tensor(cls, a, b):
        paulis = a.paulis.tensor(b.paulis)
        coeffs = np.kron(a.coeffs, b.coeffs)
        return SparsePauliOp(paulis, coeffs, copy=False)

    def _add(self, other, qargs=None):
        if qargs is None:
            qargs = getattr(other, "qargs", None)

        if not isinstance(other, SparsePauliOp):
            other = SparsePauliOp(other)

        self._op_shape._validate_add(other._op_shape, qargs)

        paulis = self.paulis._add(other.paulis, qargs=qargs)
        coeffs = np.hstack((self.coeffs, other.coeffs))
        return SparsePauliOp(paulis, coeffs, ignore_pauli_phase=True, copy=False)

    def _multiply(self, other):
        if not isinstance(other, Number):
            raise QiskitError("other is not a number")
        if other == 0:
            # Check edge case that we deleted all Paulis
            # In this case we return an identity Pauli with a zero coefficient
            paulis = PauliList.from_symplectic(
                np.zeros((1, self.num_qubits), dtype=bool),
                np.zeros((1, self.num_qubits), dtype=bool),
            )
            coeffs = np.array([0j])
            return SparsePauliOp(paulis, coeffs, ignore_pauli_phase=True, copy=False)
        # Otherwise we just update the phases
        return SparsePauliOp(
            self.paulis.copy(), other * self.coeffs, ignore_pauli_phase=True, copy=False
        )

    # ---------------------------------------------------------------------
    # Utility Methods
    # ---------------------------------------------------------------------

    def is_unitary(self, atol=None, rtol=None):
        """Return True if operator is a unitary matrix.

        Args:
            atol (float): Optional. Absolute tolerance for checking if
                          coefficients are zero (Default: 1e-8).
            rtol (float): Optional. relative tolerance for checking if
                          coefficients are zero (Default: 1e-5).

        Returns:
            bool: True if the operator is unitary, False otherwise.
        """
        # Get default atol and rtol
        if atol is None:
            atol = self.atol
        if rtol is None:
            rtol = self.rtol

        # Compose with adjoint
        val = self.compose(self.adjoint()).simplify()
        # See if the result is an identity
        return (
            val.size == 1
            and np.isclose(val.coeffs[0], 1.0, atol=atol, rtol=rtol)
            and not np.any(val.paulis.x)
            and not np.any(val.paulis.z)
        )

    def simplify(self, atol=None, rtol=None):
        """Simplify PauliList by combining duplicates and removing zeros.

        Args:
            atol (float): Optional. Absolute tolerance for checking if
                          coefficients are zero (Default: 1e-8).
            rtol (float): Optional. relative tolerance for checking if
                          coefficients are zero (Default: 1e-5).

        Returns:
            SparsePauliOp: the simplified SparsePauliOp operator.
        """
        # Get default atol and rtol
        if atol is None:
            atol = self.atol
        if rtol is None:
            rtol = self.rtol

        # Filter non-zero coefficients
        non_zero = np.logical_not(np.isclose(self.coeffs, 0, atol=atol, rtol=rtol))
        paulis_x = self.paulis.x[non_zero]
        paulis_z = self.paulis.z[non_zero]
        nz_coeffs = self.coeffs[non_zero]

        # Pack bool vectors into np.uint8 vectors by np.packbits
        array = np.packbits(paulis_x, axis=1) * 256 + np.packbits(paulis_z, axis=1)
        indexes, inverses = unordered_unique(array)

        if np.all(non_zero) and indexes.shape[0] == array.shape[0]:
            # No zero operator or duplicate operator
            return self.copy()

        coeffs = np.zeros(indexes.shape[0], dtype=complex)
        np.add.at(coeffs, inverses, nz_coeffs)
        # Delete zero coefficient rows
        is_zero = np.isclose(coeffs, 0, atol=atol, rtol=rtol)
        # Check edge case that we deleted all Paulis
        # In this case we return an identity Pauli with a zero coefficient
        if np.all(is_zero):
            x = np.zeros((1, self.num_qubits), dtype=bool)
            z = np.zeros((1, self.num_qubits), dtype=bool)
            coeffs = np.array([0j], dtype=complex)
        else:
            non_zero = np.logical_not(is_zero)
            non_zero_indexes = indexes[non_zero]
            x = paulis_x[non_zero_indexes]
            z = paulis_z[non_zero_indexes]
            coeffs = coeffs[non_zero]

        return SparsePauliOp(
            PauliList.from_symplectic(z, x), coeffs, ignore_pauli_phase=True, copy=False
        )

    def chop(self, tol=1e-14):
        """Set real and imaginary parts of the coefficients to 0 if ``< tol`` in magnitude.

        For example, the operator representing ``1+1e-17j X + 1e-17 Y`` with a tolerance larger
        than ``1e-17`` will be reduced to ``1 X`` whereas :meth:`.SparsePauliOp.simplify` would
        return ``1+1e-17j X``.

        If a both the real and imaginary part of a coefficient is 0 after chopping, the
        corresponding Pauli is removed from the operator.

        Args:
            tol (float): The absolute tolerance to check whether a real or imaginary part should
                be set to 0.

        Returns:
            SparsePauliOp: This operator with chopped coefficients.
        """
        realpart_nonzero = np.abs(self.coeffs.real) > tol
        imagpart_nonzero = np.abs(self.coeffs.imag) > tol

        remaining_indices = np.logical_or(realpart_nonzero, imagpart_nonzero)
        remaining_real = realpart_nonzero[remaining_indices]
        remaining_imag = imagpart_nonzero[remaining_indices]

        if len(remaining_real) == 0:  # if no Paulis are left
            x = np.zeros((1, self.num_qubits), dtype=bool)
            z = np.zeros((1, self.num_qubits), dtype=bool)
            coeffs = np.array([0j], dtype=complex)
        else:
            coeffs = np.zeros(np.count_nonzero(remaining_indices), dtype=complex)
            coeffs.real[remaining_real] = self.coeffs.real[realpart_nonzero]
            coeffs.imag[remaining_imag] = self.coeffs.imag[imagpart_nonzero]

            x = self.paulis.x[remaining_indices]
            z = self.paulis.z[remaining_indices]

        return SparsePauliOp(
            PauliList.from_symplectic(z, x), coeffs, ignore_pauli_phase=True, copy=False
        )

    @staticmethod
    def sum(ops):
        """Sum of SparsePauliOps.

        This is a specialized version of the builtin ``sum`` function for SparsePauliOp
        with smaller overhead.

        Args:
            ops (list[SparsePauliOp]): a list of SparsePauliOps.

        Returns:
            SparsePauliOp: the SparsePauliOp representing the sum of the input list.

        Raises:
            QiskitError: if the input list is empty.
            QiskitError: if the input list includes an object that is not SparsePauliOp.
            QiskitError: if the numbers of qubits of the objects in the input list do not match.
        """
        if len(ops) == 0:
            raise QiskitError("Input list is empty")
        if not all(isinstance(op, SparsePauliOp) for op in ops):
            raise QiskitError("Input list includes an object that is not SparsePauliOp")
        for other in ops[1:]:
            ops[0]._op_shape._validate_add(other._op_shape)

        z = np.vstack([op.paulis.z for op in ops])
        x = np.vstack([op.paulis.x for op in ops])
        phase = np.hstack([op.paulis._phase for op in ops])
        pauli_list = PauliList(BasePauli(z, x, phase))
        coeffs = np.hstack([op.coeffs for op in ops])
        return SparsePauliOp(pauli_list, coeffs, ignore_pauli_phase=True, copy=False)

    # ---------------------------------------------------------------------
    # Additional conversions
    # ---------------------------------------------------------------------

    @staticmethod
    def from_operator(obj, atol=None, rtol=None):
        """Construct from an Operator objector.

        Note that the cost of this construction is exponential as it involves
        taking inner products with every element of the N-qubit Pauli basis.

        Args:
            obj (Operator): an N-qubit operator.
            atol (float): Optional. Absolute tolerance for checking if
                          coefficients are zero (Default: 1e-8).
            rtol (float): Optional. relative tolerance for checking if
                          coefficients are zero (Default: 1e-5).

        Returns:
            SparsePauliOp: the SparsePauliOp representation of the operator.

        Raises:
            QiskitError: if the input operator is not an N-qubit operator.
        """
        # Get default atol and rtol
        if atol is None:
            atol = SparsePauliOp.atol
        if rtol is None:
            rtol = SparsePauliOp.rtol

        if not isinstance(obj, Operator):
            obj = Operator(obj)

        # Check dimension is N-qubit operator
        num_qubits = obj.num_qubits
        if num_qubits is None:
            raise QiskitError("Input Operator is not an N-qubit operator.")
        data = obj.data

        # Index of non-zero basis elements
        inds = []
        # Non-zero coefficients
        coeffs = []
        # Non-normalized basis factor
        denom = 2**num_qubits
        # Compute coefficients from basis
        basis = pauli_basis(num_qubits, pauli_list=True)
        for i, mat in enumerate(basis.matrix_iter()):
            coeff = np.trace(mat.dot(data)) / denom
            if not np.isclose(coeff, 0, atol=atol, rtol=rtol):
                inds.append(i)
                coeffs.append(coeff)
        # Get Non-zero coeff PauliList terms
        paulis = basis[inds]
        return SparsePauliOp(paulis, coeffs, copy=False)

    @staticmethod
    def from_list(obj):
        """Construct from a list of Pauli strings and coefficients.

        For example, the 5-qubit Hamiltonian

        .. math::

            H = Z_1 X_4 + 2 Y_0 Y_3

        can be constructed as

        .. code-block:: python

            # via tuples and the full Pauli string
            op = SparsePauliOp.from_list([("XIIZI", 1), ("IYIIY", 2)])

        Args:
            obj (Iterable[Tuple[str, complex]]): The list of 2-tuples specifying the Pauli terms.

        Returns:
            SparsePauliOp: The SparsePauliOp representation of the Pauli terms.

        Raises:
            QiskitError: If the list of Paulis is empty.
        """
        obj = list(obj)  # To convert zip or other iterable

        size = len(obj)  # number of Pauli terms
        if size == 0:
            raise QiskitError("Input Pauli list is empty.")

        # determine the number of qubits
        num_qubits = len(obj[0][0])

        coeffs = np.zeros(size, dtype=complex)
        labels = np.zeros(size, dtype=f"<U{num_qubits}")
        for i, item in enumerate(obj):
            labels[i] = item[0]
            coeffs[i] = item[1]

        paulis = PauliList(labels)
        return SparsePauliOp(paulis, coeffs, copy=False)

    @staticmethod
    def from_sparse_list(obj, num_qubits):
        """Construct from a list of local Pauli strings and coefficients.

        Each list element is a 3-tuple of a local Pauli string, indices where to apply it,
        and a coefficient.

        For example, the 5-qubit Hamiltonian

        .. math::

            H = Z_1 X_4 + 2 Y_0 Y_3

        can be constructed as

        .. code-block:: python

            # via triples and local Paulis with indices
            op = SparsePauliOp.from_sparse_list([("ZX", [1, 4], 1), ("YY", [0, 3], 2)], num_qubits=5)

            # equals the following construction from "dense" Paulis
            op = SparsePauliOp.from_list([("XIIZI", 1), ("IYIIY", 2)])

        Args:
            obj (Iterable[Tuple[str, List[int], complex]]): The list 3-tuples specifying the Paulis.
            num_qubits (int): The number of qubits of the operator.

        Returns:
            SparsePauliOp: The SparsePauliOp representation of the Pauli terms.

        Raises:
            QiskitError: If the list of Paulis is empty.
            QiskitError: If the number of qubits is incompatible with the indices of the Pauli terms.
        """
        obj = list(obj)  # To convert zip or other iterable

        size = len(obj)  # number of Pauli terms
        if size == 0:
            raise QiskitError("Input Pauli list is empty.")

        coeffs = np.zeros(size, dtype=complex)
        labels = np.zeros(size, dtype=f"<U{num_qubits}")

        for i, (paulis, indices, coeff) in enumerate(obj):
            # construct the full label based off the non-trivial Paulis and indices
            label = ["I"] * num_qubits
            for pauli, index in zip(paulis, indices):
                if index >= num_qubits:
                    raise QiskitError(
                        f"The number of qubits ({num_qubits}) is smaller than a required index {index}."
                    )
                label[~index] = pauli

            labels[i] = "".join(label)
            coeffs[i] = coeff

        paulis = PauliList(labels)
        return SparsePauliOp(paulis, coeffs, copy=False)

    def to_list(self, array=False):
        """Convert to a list Pauli string labels and coefficients.

        For operators with a lot of terms converting using the ``array=True``
        kwarg will be more efficient since it allocates memory for
        the full Numpy array of labels in advance.

        Args:
            array (bool): return a Numpy array if True, otherwise
                          return a list (Default: False).

        Returns:
            list or array: List of pairs (label, coeff) for rows of the PauliList.
        """
        # Dtype for a structured array with string labels and complex coeffs
        pauli_labels = self.paulis.to_labels(array=True)
        labels = np.zeros(self.size, dtype=[("labels", pauli_labels.dtype), ("coeffs", "c16")])
        labels["labels"] = pauli_labels
        labels["coeffs"] = self.coeffs
        if array:
            return labels
        return labels.tolist()

    def to_matrix(self, sparse=False):
        """Convert to a dense or sparse matrix.

        Args:
            sparse (bool): if True return a sparse CSR matrix, otherwise
                           return dense Numpy array (Default: False).

        Returns:
            array: A dense matrix if `sparse=False`.
            csr_matrix: A sparse matrix in CSR format if `sparse=True`.
        """
        mat = None
        for i in self.matrix_iter(sparse=sparse):
            if mat is None:
                mat = i
            else:
                mat += i
        return mat

    def to_operator(self):
        """Convert to a matrix Operator object"""
        return Operator(self.to_matrix())

    # ---------------------------------------------------------------------
    # Custom Iterators
    # ---------------------------------------------------------------------

    def label_iter(self):
        """Return a label representation iterator.

        This is a lazy iterator that converts each term in the SparsePauliOp
        into a tuple (label, coeff). To convert the entire table to labels
        use the :meth:`to_labels` method.

        Returns:
            LabelIterator: label iterator object for the PauliTable.
        """

        class LabelIterator(CustomIterator):
            """Label representation iteration and item access."""

            def __repr__(self):
                return f"<SparsePauliOp_label_iterator at {hex(id(self))}>"

            def __getitem__(self, key):
                coeff = self.obj.coeffs[key]
                pauli = self.obj.paulis.label_iter()[key]
                return (pauli, coeff)

        return LabelIterator(self)

    def matrix_iter(self, sparse=False):
        """Return a matrix representation iterator.

        This is a lazy iterator that converts each term in the SparsePauliOp
        into a matrix as it is used. To convert to a single matrix use the
        :meth:`to_matrix` method.

        Args:
            sparse (bool): optionally return sparse CSR matrices if True,
                           otherwise return Numpy array matrices
                           (Default: False)

        Returns:
            MatrixIterator: matrix iterator object for the PauliList.
        """

        class MatrixIterator(CustomIterator):
            """Matrix representation iteration and item access."""

            def __repr__(self):
                return f"<SparsePauliOp_matrix_iterator at {hex(id(self))}>"

            def __getitem__(self, key):
                coeff = self.obj.coeffs[key]
                mat = self.obj.paulis[key].to_matrix(sparse)
                return coeff * mat

        return MatrixIterator(self)

    def _create_graph(self, qubit_wise):
        """Transform measurement operator grouping problem into graph coloring problem

        Args:
            qubit_wise (bool): whether the commutation rule is applied to the whole operator,
                or on a per-qubit basis.

        Returns:
            retworkx.PyGraph: A class of undirected graphs
        """

        edges = self.paulis._noncommutation_graph(qubit_wise)
        graph = rx.PyGraph()
        graph.add_nodes_from(range(self.size))
        graph.add_edges_from_no_data(edges)
        return graph

    def group_commuting(self, qubit_wise=False):
        """Partition a SparsePauliOp into sets of commuting Pauli strings.

        Args:
            qubit_wise (bool): whether the commutation rule is applied to the whole operator,
                or on a per-qubit basis.  For example:

                .. code-block:: python

                    >>> op = SparsePauliOp.from_list([("XX", 2), ("YY", 1), ("IZ",2j), ("ZZ",1j)])
                    >>> op.group_commuting()
                    [SparsePauliOp(["IZ", "ZZ"], coeffs=[0.+2.j, 0.+1j]),
                     SparsePauliOp(["XX", "YY"], coeffs=[2.+0.j, 1.+0.j])]
                    >>> op.group_commuting(qubit_wise=True)
                    [SparsePauliOp(['XX'], coeffs=[2.+0.j]),
                     SparsePauliOp(['YY'], coeffs=[1.+0.j]),
                     SparsePauliOp(['IZ', 'ZZ'], coeffs=[0.+2.j, 0.+1.j])]

        Returns:
            List[SparsePauliOp]: List of SparsePauliOp where each SparsePauliOp contains
                commuting Pauli operators.
        """

        graph = self._create_graph(qubit_wise)
        # Keys in coloring_dict are nodes, values are colors
        coloring_dict = rx.graph_greedy_color(graph)
        groups = defaultdict(list)
        for idx, color in coloring_dict.items():
            groups[color].append(idx)
        return [self[group] for group in groups.values()]


# ibm-quantum-spring-challenge-2022/exercise4/04.quantum_chemistry.ipynb

from qiskit_nature.drivers import Molecule
from qiskit_nature.drivers.second_quantization import ElectronicStructureMoleculeDriver, ElectronicStructureDriverType

molecule = Molecule(
    # coordinates are given in Angstrom
    geometry=[
        ["O", [0.0, 0.0, 0.0]],
        ["H", [0.758602, 0.0, 0.504284]],
        ["H", [0.758602, 0.0, -0.504284]]
    ],
    multiplicity=1,  # = 2*spin + 1
    charge=0,
)

driver = ElectronicStructureMoleculeDriver(
    molecule=molecule,
    basis="sto3g",
    driver_type=ElectronicStructureDriverType.PYSCF,
)

properties = driver.run()



from qiskit_nature.transformers.second_quantization.electronic import ActiveSpaceTransformer

# Check the occupation of the spin orbitals
PN_property = properties.get_property("ParticleNumber")
print(PN_property)

# Define the active space around the Fermi level 
# (selected automatically around the HOMO and LUMO, ordered by energy)
transformer = ActiveSpaceTransformer(
    num_electrons=2, #how many electrons we have in our active space
    num_molecular_orbitals=2, #how many orbitals we have in our active space
)


# We can hand-pick the MOs to be included in the AS
# (in case they are not exactly around the Fermi level)

# transformer = ActiveSpaceTransformer(
#     num_electrons=2, #how amny electrons we have in our active space
#     num_molecular_orbitals=2, #how many orbitals we have in our active space
#     active_orbitals=[4,7], #We are specifically choosing MO number 4 (occupied with two electrons) and 7 (empty)
# )

from qiskit_nature.problems.second_quantization.electronic import ElectronicStructureProblem

problem = ElectronicStructureProblem(driver, [transformer])

second_q_ops = problem.second_q_ops() # this calls driver.run() internally

hamiltonian = second_q_ops[0]

print(hamiltonian)

from qiskit_nature.mappers.second_quantization import ParityMapper, BravyiKitaevMapper, JordanWignerMapper
from qiskit_nature.converters.second_quantization import QubitConverter

# Setup the mapper and qubit converter
mapper_type = 'JordanWignerMapper'

if mapper_type == 'ParityMapper':
    mapper = ParityMapper()
elif mapper_type == 'JordanWignerMapper':
    mapper = JordanWignerMapper()
elif mapper_type == 'BravyiKitaevMapper':
    mapper = BravyiKitaevMapper()

converter = QubitConverter(mapper)

qubit_op = converter.convert(hamiltonian)
print(qubit_op)

print("No-grouping (Active Space)")
print("\n", f"Number of Measurement is {len(qubit_op)}")

print("Qubit-wise grouping")
measure_group = PauliList(qubit_op.primitive.paulis).group_commuting(qubit_wise=True)
for group in measure_group:
    print(group)
    
print("\n", f"Number of Measurement is {len(measure_group)}")

print("General grouping")
measure_group = PauliList(qubit_op.primitive.paulis).group_commuting()
for group in measure_group:
    print(group)
    
print("\n", f"Number of Measurement is {len(measure_group)}")

problem = ElectronicStructureProblem(driver, [])

second_q_ops = problem.second_q_ops() # this calls driver.run() internally

hamiltonian = second_q_ops[0]

# print(hamiltonian)
converter = QubitConverter(mapper)

qubit_op = converter.convert(hamiltonian)
print(qubit_op)

print("No-grouping (Active Space)")
print("\n", f"Number of Measurement is {len(qubit_op)}")

print("Qubit-wise grouping")
measure_group = PauliList(qubit_op.primitive.paulis).group_commuting(qubit_wise=True)
print("\n", f"Number of Measurement is {len(measure_group)}\n")

for group in measure_group:
    print(group)
    


print("General grouping")
measure_group = PauliList(qubit_op.primitive.paulis).group_commuting()
print("\n", f"Number of Measurement is {len(measure_group)}\n")
for group in measure_group:
    print(group)
    




