fix: 🐛 Making number of qubits obligatory

benches/_builders.py

@@ -167,6 +167,7 @@ def make_random_problem(
         majoranas=gen_majoranas,
         gen_coeffs=gen_coeffs,
         param_inds=param_inds,
+        system_size=num_modes,
         parameters=parameters,
         initial_state=[],
     )
@@ -293,7 +294,12 @@ def build_hubbard_problem(
     fermi_gates = [ExpGate(term) for term in _hubbard_fermion_terms(config)]
     parameters = [config.trotter_dt] * len(fermi_gates)
     occupied = _neel_occupied_modes(config.num_sites, config.neel_start_spin)
-    circuit = Circuit(gates=fermi_gates, parameters=parameters, initial_state=occupied)
+    circuit = Circuit(
+        gates=fermi_gates,
+        parameters=parameters,
+        initial_state=occupied,
+        system_size=config.num_qubits,
+    )
 
     observable = FermiOperator(
         terms=[
@@ -525,6 +531,7 @@ def build_kicked_ising_problem(
         gates=tuple(gate for gate, _ in gate_angles),
         parameters=tuple(angle for _, angle in gate_angles),
         initial_state=[],
+        system_size=config.num_qubits,
     )
 
     obs_str = (

benches/third_party/majorana_prop/monoprop_hubbard1d_benchmark.py

@@ -157,6 +157,7 @@ def main():
         gates=trotter_gates,
         parameters=trotter_parameters,
         initial_state=intial_state,
+        system_size=num_qubits,
     )
 
     simulator = MajoranaPropagator(

src/monoprop/circuit.py

@@ -57,6 +57,7 @@
 from .conversion_utils import _extend_pauli_string, _pauli_to_majorana
 from .majorana import MajoranaOperator
 from .pauli import Pauli, PauliOperator
+from .utils import _validate_system_size
 
 if TYPE_CHECKING:
     from collections.abc import Iterator, Sequence
@@ -223,6 +224,21 @@ def __repr__(self) -> str:
         """Return a string representation such as ``ExpGate(<generator>, index=0)``."""
         return f"{self.__class__.__name__}({self.generator}, index={self.index})"
 
+    @property
+    def system_size(self) -> int:
+        """Number of modes/qubits the generator acts on.
+
+        Reads ``num_modes`` off a [MajoranaOperator][monoprop.majorana.MajoranaOperator]
+        generator, or ``num_qubits`` off a [PauliOperator][monoprop.pauli.PauliOperator]
+        generator.
+
+        Returns:
+            Number of modes/qubits
+        """
+        if isinstance(self.generator, PauliOperator):
+            return self.generator.num_qubits
+        return self.generator.num_modes
+
 
 class Circuit:
     """A variational circuit: an ordered sequence of exponential gates, angles, and a state.
@@ -255,33 +271,42 @@ class Circuit:
         gates: The ordered exponential gates.
         parameters: The angle values, or empty for an unbound circuit.
         initial_state: The reference state (occupied mode / qubit indices).
+        system_size: System width (number of fermionic modes / qubits).
         family: The gate family -- ``"pauli"``, ``"majorana"``, or ``"empty"`` -- computed at
             construction; the propagators dispatch on it.
     """
 
-    def __init__(
+    def __init__(  # noqa:C901
         self,
-        gates: Sequence[ExpGate] = (),
-        parameters: Sequence[float] = (),
+        gates: Sequence[ExpGate],
+        system_size: int,
         initial_state: Sequence[int] = (),
+        parameters: Sequence[float] = (),
     ) -> None:
         """Build the circuit, dropping identity gates and validating family/mapping/params.
 
         Args:
             gates: The ordered exponential gates.
-            parameters: The angle values, or empty for an unbound circuit.
             initial_state: The reference state (occupied mode / qubit indices).
+            system_size: Number of modes/qubits for the circuit.
+            parameters: The angle values, or empty for an unbound circuit.
 
         Raises:
             ValueError: On duplicate initial-state indices, a bad parameter mapping, or a
-                bound circuit whose parameter count does not match [n_parameters][].
+                bound circuit whose parameter count does not match [n_parameters][]; also if
+                a gate's operator width differs from ``system_size``.
             TypeError: On a non-[ExpGate][] gate or a mix of qubit and Majorana gate families.
         """
         gates = tuple(gates)
         parameters = tuple(float(v) for v in parameters)
         initial_state = tuple(int(i) for i in initial_state)
+        system_size = _validate_system_size(system_size, argument_name="system_size")
         if len(set(initial_state)) != len(initial_state):
             raise ValueError("Duplicate indices in initial state")
+        if any(i < 0 or i >= system_size for i in initial_state):
+            raise ValueError(
+                f"initial_state entries must be in 0..{system_size - 1}; got {list(initial_state)}."
+            )
 
         # Validate gate types up front: the identity-drop below reads gate.index/.generator,
         # so a non-ExpGate gate must be rejected with a clear TypeError first rather than crashing
@@ -291,6 +316,11 @@ def __init__(
                 raise TypeError(
                     f"Circuit gates must be ExpGate; got {type(gate).__name__}."
                 )
+            gate_size = gate.system_size
+            if gate_size != system_size:
+                raise ValueError(
+                    f"Gate generator width {gate_size} does not match circuit system_size={system_size}."
+                )
 
         def _is_identity_gate(gate: ExpGate) -> bool:
             return all(coeff == 0 for coeff in gate.generator.terms.values())
@@ -310,6 +340,7 @@ def _is_identity_gate(gate: ExpGate) -> bool:
         self.gates = gates
         self.parameters = parameters
         self.initial_state = initial_state
+        self.system_size = system_size
         #: The gate family, computed from the (validated) gates; the propagators dispatch on it.
         self.family = self._resolve_family(gates)
 
@@ -328,6 +359,7 @@ def __eq__(self, other: object) -> bool:
             self.gates == other.gates
             and self.parameters == other.parameters
             and self.initial_state == other.initial_state
+            and self.system_size == other.system_size
         )
 
     __hash__ = None  # type: ignore[assignment]  # value-equal but not hashable (mutable gates)
@@ -336,7 +368,8 @@ def __repr__(self) -> str:
         """Return a string representation listing the gates, parameters, and initial state."""
         return (
             f"{self.__class__.__name__}(gates={self.gates!r}, "
-            f"parameters={self.parameters!r}, initial_state={self.initial_state!r})"
+            f"parameters={self.parameters!r}, initial_state={self.initial_state!r}, "
+            f"system_size={self.system_size!r})"
         )
 
     @staticmethod
@@ -420,6 +453,11 @@ def __add__(self, other: Circuit) -> Circuit:
             raise ValueError(
                 "Cannot concatenate circuits with different initial states."
             )
+        if self.system_size != other.system_size:
+            raise ValueError(
+                f"Cannot concatenate circuits with different system_size: "
+                f"{self.system_size} != {other.system_size}."
+            )
         offset = self.n_parameters
         # Preserve each gate's _structural flag: a dense (wire-format) gate carries structural
         # coefficients that must not be antihermitian-normalized again.
@@ -435,6 +473,7 @@ def __add__(self, other: Circuit) -> Circuit:
             gates=left + right,
             parameters=tuple(self.parameters) + tuple(other.parameters),
             initial_state=self.initial_state or other.initial_state,
+            system_size=self.system_size,
         )
 
     @classmethod
@@ -443,6 +482,7 @@ def from_dense_arrays(
         majoranas: Sequence[Sequence[int]],
         gen_coeffs: Sequence[float],
         param_inds: Sequence[int],
+        system_size: int,
         parameters: Sequence[float] = (),
         initial_state: Sequence[int] = (),
     ) -> Circuit:
@@ -460,13 +500,15 @@ def from_dense_arrays(
             gen_coeffs: Generator coefficient per monomial.
             param_inds: Variational-angle index per monomial (contiguous runs group into
                 gates).
+            system_size: Number of fermionic modes in the system.
             parameters: Optional angle values.
             initial_state: Optional reference state (occupied mode indices).
 
         Returns:
             A [Circuit][] carrying the grouped gates, angle values, and initial state.
         """
         indices = [int(p) for p in param_inds]
+        system_size = _validate_system_size(system_size, argument_name="system_size")
         gates: list[ExpGate] = []
         current_index: int | None = None
         current_majoranas: list[tuple[int, ...]] = []
@@ -478,7 +520,7 @@ def _flush() -> None:
             gates.append(
                 ExpGate._structural_gate(
                     MajoranaOperator._from_terms(
-                        current_majoranas, current_coeffs, num_modes=None
+                        current_majoranas, current_coeffs, num_modes=system_size
                     ),
                     index=current_index,
                 )
@@ -490,7 +532,16 @@ def _flush() -> None:
                 current_majoranas = []
                 current_coeffs = []
             current_index = pidx
-            current_majoranas.append(tuple(int(i) for i in maj))
+            majorana = tuple(int(i) for i in maj)
+            if any(i < 0 for i in majorana):
+                raise ValueError(
+                    f"Majorana indices must be non-negative; got {majorana}."
+                )
+            if majorana and max(majorana) >= 2 * system_size:
+                raise ValueError(
+                    f"Majorana term {majorana} acts on an index >= 2*system_size={2 * system_size}."
+                )
+            current_majoranas.append(majorana)
             current_coeffs.append(complex(float(coeff)))
         if current_majoranas:
             _flush()
@@ -499,6 +550,7 @@ def _flush() -> None:
             gates=tuple(gates),
             parameters=tuple(float(p) for p in parameters),
             initial_state=tuple(int(i) for i in initial_state),
+            system_size=system_size,
         )
 
 
@@ -632,9 +684,7 @@ def _validate_commuting_majorana_generator(generator: MajoranaOperator) -> None:
             )
 
 
-def _gate_layers(
-    gate: ExpGate, num_qubits: int | None
-) -> list[tuple[tuple[int, ...], float]]:
+def _gate_layers(gate: ExpGate, num_qubits: int) -> list[tuple[tuple[int, ...], float]]:
     r"""Expand one gate into ``(majorana, gen_coeff)`` layers, in application order.
 
     A ``"pauli"``-family [ExpGate][] places each [Pauli][monoprop.pauli.Pauli] term on
@@ -649,8 +699,6 @@ def _gate_layers(
     # ``isinstance`` narrows the fall-through arm to MajoranaOperator).
     generator = gate.generator
     if isinstance(generator, PauliOperator):
-        if num_qubits is None:
-            raise ValueError("num_qubits is required to expand a Pauli gate.")
         layers: list[tuple[tuple[int, ...], float]] = []
         for pauli, coeff in generator.terms.items():
             extended = _extend_pauli_string(pauli.string, pauli.qubits, num_qubits)
@@ -673,7 +721,7 @@ def _gate_layers(
 def expand_monomials(
     gates: Sequence[ExpGate],
     mapping: Sequence[int],
-    num_qubits: int | None = None,
+    num_qubits: int,
 ) -> tuple[list[tuple[int, ...]], list[float], list[int], list[int]]:
     """Flatten gates + an already-resolved per-gate mapping into per-monomial arrays.
 
@@ -689,6 +737,7 @@ def expand_monomials(
         the authoring gate monomial ``i`` came from, so the engine can recover gate
         boundaries; monomials from a multi-term gate share one gate index.
     """
+    num_qubits = _validate_system_size(num_qubits, argument_name="num_qubits")
     majoranas: list[tuple[int, ...]] = []
     gen_coeffs: list[float] = []
     per_monomial: list[int] = []

src/monoprop/fermi.py

@@ -22,6 +22,7 @@
 
 from .conversion_utils import _n_product
 from .majorana import MajoranaOperator
+from .utils import _validate_system_size
 
 if TYPE_CHECKING:
     from collections.abc import Sequence
@@ -109,39 +110,43 @@ def __init__(
         self,
         terms: Sequence[FermiString] | Sequence[Sequence[tuple[int, str]]],
         coefficients: Sequence[complex],
-        num_modes: int | None = None,
+        num_modes: int,
     ) -> None:
         """Initialize the fermi operator.
 
         Args:
             terms: List of FermiString objects representing the operator.
             coefficients: List of coefficients corresponding to the terms.
-            num_modes: Optional number of modes. If not provided, it will be inferred from the terms.
+            num_modes: Number of modes in the system.
 
         Raises:
-            ValueError: If any index is out of bounds or if there are duplicate indices.
+            TypeError: If ``num_modes`` is not an integer.
+            ValueError: If ``num_modes`` is negative or any index is out of bounds.
         """
         self.terms = [
             t if isinstance(t, FermiString) else FermiString(t) for t in terms
         ]
         self.coefficients = list(coefficients)
-        self.num_modes = (
-            num_modes
-            if num_modes is not None
-            else max((idx for f in self.terms for idx, _ in f.expression)) + 1
-        )
+        self.num_modes = _validate_system_size(num_modes, argument_name="num_modes")
+        for term in self.terms:
+            for idx, _ in term.expression:
+                if idx >= self.num_modes:
+                    raise ValueError(
+                        "Fermi term index out of bounds: "
+                        f"{idx} >= num_modes={self.num_modes}."
+                    )
 
     @classmethod
     def from_dict(
-        cls, terms_dict: dict[tuple[tuple[int, str], ...], complex]
+        cls, terms_dict: dict[tuple[tuple[int, str], ...], complex], num_modes: int
     ) -> FermiOperator:
         """Construct a FermiOperator from a dictionary."""
         terms = []
         coefficients = []
         for key, value in terms_dict.items():
             terms.append(FermiString(key))
             coefficients.append(value)
-        return cls(terms=terms, coefficients=coefficients)
+        return cls(terms=terms, coefficients=coefficients, num_modes=num_modes)
 
     def __len__(self) -> int:
         """Number of terms in the operator."""

src/monoprop/integral_conversion.py

@@ -104,4 +104,5 @@ def integrals_to_fermion(
             continue
         terms[ind] += coeff
 
-    return FermiOperator.from_dict(terms)
+    # 2* because we have alpha and beta orbitals
+    return FermiOperator.from_dict(terms, 2 * hamiltonian[1].shape[1])

src/monoprop/majorana.py

@@ -21,6 +21,8 @@
 
 import numpy as np
 
+from .utils import _validate_system_size
+
 if TYPE_CHECKING:
     from collections.abc import Mapping, Sequence
 
@@ -99,30 +101,42 @@ def __init__(
                 also authored as a [MajoranaOperator][] (wrapped in
                 [ExpGate][monoprop.circuit.ExpGate]) -- bare [Majorana][] terms are not accepted
                 by ``ExpGate``, since the operator is what carries the mode count.
+
+        Raises:
+            TypeError: If ``num_modes`` is not an integer.
+            ValueError: If ``num_modes`` is negative or a term index is out of range.
         """
         # Route raw index tuples through Majorana so they get the same non-negative/distinct
         # validation a Majorana key already carries (a bare tuple would otherwise slip past it).
         majoranas = [
             (key if isinstance(key, Majorana) else Majorana(*key)).indices
             for key in terms
         ]
-        self.num_modes = num_modes
+        self.num_modes = _validate_system_size(num_modes, argument_name="num_modes")
+        for majorana in majoranas:
+            # majoranas are sorted in here
+            if majorana and majorana[-1] >= 2 * self.num_modes:
+                raise ValueError(
+                    f"Majorana term {majorana} acts on an index >= num_modes={self.num_modes}."
+                )
         self.terms = self._accumulate(majoranas, list(terms.values()))
 
     @classmethod
     def _from_terms(
         cls,
         majoranas: Sequence[Sequence[int]],
         coefficients: Sequence[complex],
-        num_modes: int | None = None,
+        num_modes: int,
     ) -> MajoranaOperator:
         """Build from parallel ``majoranas``/``coefficients`` lists (internal).
 
         Unlike the dict constructor this accepts colliding monomials and sums them, which the
         Jordan-Wigner and fermionic conversions ([get_majorana_operator][]) rely on.
         """
         obj = cls.__new__(cls)
-        obj.num_modes = num_modes
+        obj.num_modes = _validate_system_size(num_modes, argument_name="num_modes")
+        # _accumulate is taking care of the sorting/removing repeating terms
+        # so no need to check in here
         obj.terms = cls._accumulate(majoranas, coefficients)
         return obj