Implement SimulatorBackend and Statevec

docs/source/matrix.rst

@@ -10,3 +10,5 @@ Functions
 ---------
 
 .. autofunction:: graphix_zx.matrix.is_unitary
+
+.. autofunction:: graphix_zx.matrix.is_hermitian

docs/source/references.rst

@@ -5,6 +5,8 @@ Module reference
     :maxdepth: 2
 
     common
+    simulator_backend
+    statevec
     euler
     matrix
     graphstate

docs/source/simulator_backend.rst

@@ -0,0 +1,7 @@
+Simulator Backend
+=================
+
+.. automodule:: graphix_zx.simulator_backend
+   :members:
+   :undoc-members:
+   :show-inheritance:

docs/source/statevec.rst

@@ -0,0 +1,7 @@
+State Vector
+============
+
+.. automodule:: graphix_zx.statevec
+   :members:
+   :undoc-members:
+   :show-inheritance:

graphix_zx/common.py

@@ -70,6 +70,16 @@ def angle(self) -> float:
         """Return the measurement angle."""
         raise NotImplementedError
 
+    @abc.abstractmethod
+    def flip(self) -> MeasBasis:
+        """Flip the measurement basis."""
+        raise NotImplementedError
+
+    @abc.abstractmethod
+    def conjugate(self) -> MeasBasis:
+        """Return the conjugate of the measurement basis."""
+        raise NotImplementedError
+
     @abc.abstractmethod
     def vector(self) -> NDArray[np.complex128]:
         """Return the measurement basis vector."""
@@ -108,16 +118,17 @@ def angle(self) -> float:
         return self.__angle
 
     @typing_extensions.override
-    def vector(self) -> NDArray[np.complex128]:
-        r"""Return the measurement basis vector.
+    def flip(self) -> PlannerMeasBasis:
+        """Flip the measurement basis.
 
         Returns
         -------
-        `numpy.typing.NDArray`\[`numpy.complex128`\]
-            measurement basis vector
+        `PlannerMeasBasis`
+            flipped PlannerMeasBasis
         """
-        return meas_basis(self.plane, self.angle)
+        return PlannerMeasBasis(self.plane, self.angle + np.pi)
 
+    @typing_extensions.override
     def conjugate(self) -> PlannerMeasBasis:
         """Return the conjugate of the PlannerMeasBasis object.
 
@@ -142,6 +153,17 @@ def conjugate(self) -> PlannerMeasBasis:
             return PlannerMeasBasis(Plane.XZ, self.angle)
         typing_extensions.assert_never(self.plane)
 
+    @typing_extensions.override
+    def vector(self) -> NDArray[np.complex128]:
+        r"""Return the measurement basis vector.
+
+        Returns
+        -------
+        `numpy.typing.NDArray`\[`numpy.complex128`\]
+            measurement basis vector
+        """
+        return meas_basis(self.plane, self.angle)
+
 
 class AxisMeasBasis(MeasBasis):
     """Class to represent an axis measurement basis.
@@ -212,6 +234,30 @@ def angle(self) -> float:
             angle = 0 if self.sign == Sign.PLUS else np.pi
         return angle
 
+    @typing_extensions.override
+    def flip(self) -> AxisMeasBasis:
+        """Flip the measurement basis.
+
+        Returns
+        -------
+        `AxisMeasBasis`
+            flipped AxisMeasBasis
+        """
+        return AxisMeasBasis(self.axis, Sign.MINUS if self.sign == Sign.PLUS else Sign.PLUS)
+
+    @typing_extensions.override
+    def conjugate(self) -> AxisMeasBasis:
+        """Return the conjugate of the AxisMeasBasis object.
+
+        Returns
+        -------
+        `AxisMeasBasis`
+            conjugate AxisMeasBasis
+        """
+        if self.axis == Axis.Y:
+            return AxisMeasBasis(Axis.Y, Sign.MINUS if self.sign == Sign.PLUS else Sign.PLUS)
+        return AxisMeasBasis(self.axis, self.sign)
+
     @typing_extensions.override
     def vector(self) -> NDArray[np.complex128]:
         r"""Return the measurement basis vector.

graphix_zx/matrix.py

@@ -3,6 +3,7 @@
 This module provides:
 
 - `is_unitary`: check if a matrix is unitary.
+- `is_hermitian`: check if a matrix is Hermitian.
 """
 
 from __future__ import annotations
@@ -33,3 +34,19 @@ def is_unitary(mat: NDArray[T]) -> bool:
     if mat.shape[0] != mat.shape[1]:
         return False
     return np.allclose(np.eye(mat.shape[0]), mat @ mat.T.conj())
+
+
+def is_hermitian(mat: NDArray[T]) -> bool:
+    r"""Check if a matrix is Hermitian.
+
+    Parameters
+    ----------
+    mat : `numpy.typing.NDArray`\[T\]
+        matrix to check
+
+    Returns
+    -------
+    `bool`
+        `True` if Hermitian, `False` otherwise
+    """
+    return np.allclose(mat, mat.T.conj())

graphix_zx/simulator_backend.py

@@ -0,0 +1,205 @@
+"""The base class for simulator backends.
+
+This module provides:
+
+- `QubitIndexManager`: Manages the mapping of external qubit indices to internal indices
+- `BaseSimulatorBackend`: Abstract base class for simulator backends.
+
+"""
+
+from __future__ import annotations
+
+import abc
+import itertools
+import typing
+from abc import ABC
+from typing import TYPE_CHECKING
+
+if TYPE_CHECKING:
+    from collections.abc import Sequence
+
+    import numpy as np
+    from numpy.typing import NDArray
+
+    from graphix_zx.common import MeasBasis
+
+
+class QubitIndexManager:
+    """Manages the mapping of external qubit indices to internal indices."""
+
+    def __init__(self, num_qubits: int) -> None:
+        """Initialize the QubitIndexManager with a list of initial indices."""
+        self.__indices = list(range(num_qubits))
+
+    @property
+    def num_qubits(self) -> int:
+        """Get the number of qubits managed by this manager.
+
+        Returns
+        -------
+        `int`
+            The number of qubits.
+        """
+        return len(self.__indices)
+
+    def add_qubits(self, num_qubits: int) -> None:
+        """Add a specified number of qubits to the index manager.
+
+        Parameters
+        ----------
+        num_qubits : `int`
+            The number of qubits to add.
+        """
+        current_max = max(self.__indices, default=-1)
+        self.__indices.extend(range(current_max + 1, current_max + 1 + num_qubits))
+
+    def remove_qubit(self, qubit: int) -> None:
+        r"""Remove specified qubit from the index manager.
+
+        Parameters
+        ----------
+        qubit : `int`
+            The qubit to remove.
+        """
+        self.__indices = [q if q < qubit else q - 1 for q in self.__indices if q != qubit]
+
+    def match(self, order: Sequence[int]) -> bool:
+        r"""Check if the current indices match the given order.
+
+        Parameters
+        ----------
+        order : `collections.abc.Sequence`\[`int`\]
+            A sequence of indices to compare against the current indices.
+
+        Returns
+        -------
+        `bool`
+            True if the current indices match the given order, False otherwise.
+        """
+        return all(lhs == rhs for lhs, rhs in itertools.zip_longest(self.__indices, order, fillvalue=None))
+
+    def reorder(self, permutation: Sequence[int]) -> None:
+        r"""Reorder the indices based on a given permutation.
+
+        if permutation is [2, 0, 1], then
+        # [q0, q1, q2] -> [q1, q2, q0]
+
+        Parameters
+        ----------
+        permutation : `collections.abc.Sequence`\[`int`\]
+            A sequence of indices that defines the new order of the indices.
+
+        Raises
+        ------
+        ValueError
+            If the length of the permutation does not match the number of indices.
+        """
+        if len(permutation) != len(self.__indices):
+            msg = "Permutation length must match the number of indices."
+            raise ValueError(msg)
+        self.__indices = [self.__indices[i] for i in permutation]
+
+    def inverse_permutation(self) -> list[int]:
+        r"""Get the permutation that would recover the original order of indices.
+
+        Returns
+        -------
+        `list`\[`int`\]
+            A sequence of indices that maps the current order back to the original order.
+        """
+        inverse_perm = [0] * len(self.__indices)
+        for i, index in enumerate(self.__indices):
+            inverse_perm[index] = i
+        return inverse_perm
+
+    @typing.overload
+    def external_to_internal(self, external_qubits: int) -> int: ...
+
+    @typing.overload
+    def external_to_internal(self, external_qubits: Sequence[int]) -> tuple[int, ...]: ...
+
+    def external_to_internal(self, external_qubits: int | Sequence[int]) -> int | tuple[int, ...]:
+        r"""Convert external qubit indices to internal indices.
+
+        Parameters
+        ----------
+        external_qubits : `int` | `collections.abc.Sequence`\[`int`\]
+            A sequence of external qubit indices.
+
+        Returns
+        -------
+        `int` | `tuple`\[`int`, ...\]
+            A list of internal qubit indices corresponding to the external ones.
+        """
+        if isinstance(external_qubits, int):
+            return self.__indices[external_qubits]
+        return tuple(self.__indices[q] for q in external_qubits)
+
+    @typing.overload
+    def internal_to_external(self, internal_qubits: int) -> int: ...
+
+    @typing.overload
+    def internal_to_external(self, internal_qubits: Sequence[int]) -> tuple[int, ...]: ...
+
+    def internal_to_external(self, internal_qubits: int | Sequence[int]) -> int | tuple[int, ...]:
+        r"""Convert internal qubit indices to external indices.
+
+        Parameters
+        ----------
+        internal_qubits : `int` | `collections.abc.Sequence`\[`int`\]
+            A sequence of internal qubit indices.
+
+        Returns
+        -------
+        `int` | `tuple`\[`int`, ...\]
+            A list of external qubit indices corresponding to the internal ones.
+        """
+        inverse_perm = self.inverse_permutation()
+        if isinstance(internal_qubits, int):
+            return inverse_perm[internal_qubits]
+        return tuple(inverse_perm[q] for q in internal_qubits)
+
+
+# backend for all simulator backends
+class BaseSimulatorBackend(ABC):
+    """Base class for simulator backends."""
+
+    @property
+    @abc.abstractmethod
+    def num_qubits(self) -> int:
+        """Get the number of qubits in the state.
+
+        Returns
+        -------
+        `int`
+            The number of qubits in the state.
+        """
+        raise NotImplementedError
+
+    @abc.abstractmethod
+    def evolve(self, operator: NDArray[np.complex128], qubits: int | Sequence[int]) -> None:
+        r"""Evolve the state by applying an operator to a subset of qubits.
+
+        Parameters
+        ----------
+        operator : `numpy.typing.NDArray`\[`numpy.complex128`\]
+            The operator to apply.
+        qubits : `int` | `collections.abc.Sequence`\[`int`\]
+            The qubits to apply the operator to.
+        """
+        raise NotImplementedError
+
+    @abc.abstractmethod
+    def measure(self, qubit: int, meas_basis: MeasBasis, result: int) -> None:
+        """Measure a qubit in a given measurement basis.
+
+        Parameters
+        ----------
+        qubit : `int`
+            The qubit to measure.
+        meas_basis : `MeasBasis`
+            The measurement basis to use.
+        result : `int`
+            The measurement result.
+        """
+        raise NotImplementedError