diff --git a/docs/source/library_map.rst b/docs/source/library_map.rst index 72c2056a7..e0b5af70f 100644 --- a/docs/source/library_map.rst +++ b/docs/source/library_map.rst @@ -133,6 +133,7 @@ Distances across packages Exact distances of binary codes use an optimized `Brouwer--Zimmermann search `_, and nonbinary codes use exhaustive enumeration. Binary codes with sparse parity checks can instead use ``method="connected_support"``, a `search over supports that are connected in the Tanner graph `_. All exact methods can take exponential time. +A code that encodes nothing, or a circuit with no undetected logical fault, has no distance, and distance methods then return :data:`~qldpc.codes.distance.core.NO_DISTANCE`, a float ``nan`` (see :class:`~qldpc.codes.distance.core.NoDistance`). ``get_distance_exact(cutoff=c)`` stops once it finds a nontrivial word of weight at most ``c``. A result of at most ``c`` is therefore only an upper bound, while a larger result, or any result with ``cutoff <= 1`` (the default is 1), is exact. Exact searches choose their routes by cost models whose settings and constants :doc:`exact_search_tuning` maps and documents. diff --git a/src/qldpc/circuits/certification/fault_distance.py b/src/qldpc/circuits/certification/fault_distance.py index 2e55043aa..8f4d800b4 100644 --- a/src/qldpc/circuits/certification/fault_distance.py +++ b/src/qldpc/circuits/certification/fault_distance.py @@ -4,12 +4,14 @@ from __future__ import annotations +from typing import Literal + import galois import numpy as np import stim from qldpc import codes, decoders -from qldpc.codes.distance import DistanceMethod +from qldpc.codes.distance import DistanceMethod, NoDistance from qldpc.objects import Pauli #################################################################################################### @@ -18,7 +20,7 @@ def get_fault_distance( circuit: stim.Circuit, *, method: DistanceMethod = "brouwer_zimmermann" -) -> int | float: +) -> int | Literal[NoDistance.NAN]: r"""Compute the exact fault distance of a noisy circuit. Consider a noisy Stim circuit annotated with detectors and logical observables. Calling @@ -71,8 +73,8 @@ def get_fault_distance( detectors (see :meth:`qldpc.codes.CSSCode.get_distance_exact`). Returns: - The minimum number of DEM error mechanisms causing an undetected observable flip, or np.nan - if no such fault exists. + The minimum number of DEM error mechanisms causing an undetected observable flip, or + ``NO_DISTANCE`` (see :class:`~qldpc.codes.distance.NoDistance`) if no such fault exists. """ dem_arrays = decoders.DetectorErrorModelArrays(circuit) detectors = galois.GF2(dem_arrays.detector_flip_matrix.toarray()) diff --git a/src/qldpc/circuits/certification/fault_distance_test.py b/src/qldpc/circuits/certification/fault_distance_test.py index d1e22b6cb..740c2881b 100644 --- a/src/qldpc/circuits/certification/fault_distance_test.py +++ b/src/qldpc/circuits/certification/fault_distance_test.py @@ -4,11 +4,10 @@ from __future__ import annotations -import numpy as np import pytest import stim -from qldpc import circuits +from qldpc import circuits, codes from qldpc.codes.distance import DistanceMethod @@ -32,6 +31,6 @@ def test_fault_distance(measurement: int, expected_distance: int, method: Distan def test_fault_distance_without_logical_faults() -> None: - """No undetected observable flip has the code distance convention np.nan.""" + """Without an undetected observable flip, a circuit has no distance, and gives NO_DISTANCE.""" circuit = stim.Circuit("R 0\nX_ERROR(0.01) 0\nM 0\nDETECTOR rec[-1]") - assert np.isnan(circuits.get_fault_distance(circuit)) + assert circuits.get_fault_distance(circuit) is codes.NO_DISTANCE diff --git a/src/qldpc/codes/__init__.py b/src/qldpc/codes/__init__.py index 128e19f8a..2ac6bc77b 100644 --- a/src/qldpc/codes/__init__.py +++ b/src/qldpc/codes/__init__.py @@ -27,11 +27,13 @@ QuditCode, ) from .distance import ( + NO_DISTANCE, CSSLogicalCounter, CSSLogicalCounts, DistanceBackend, DistanceMethod, LogicalCountMethod, + NoDistance, count_logical_ops_css, get_distance_classical, get_distance_quantum, @@ -74,6 +76,7 @@ ) __all__ = [ + "NO_DISTANCE", "AbstractCode", "BBCode", "BCHCode", @@ -106,6 +109,7 @@ "LogicalCountMethod", "MBCode", "ManyHypercubeCode", + "NoDistance", "PairPartitionCode", "QTCode", "QuantumGolayCode", diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index 9f6179033..4db15ce94 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -47,8 +47,10 @@ from ._code_capacity import CodeCapacityDecoder, get_code_capacity_decoder from .distance import ( + NO_DISTANCE, DistanceBackend, DistanceMethod, + NoDistance, get_distance_classical, get_distance_quantum, ) @@ -108,8 +110,9 @@ class AbstractCode(abc.ABC): _matrix: galois.FieldArray _field: type[galois.FieldArray] _dimension: int | None = None - _distance: int | float | None = None - _distance_bound: int | float | None = None # best known upper bound on distance + _distance: int | Literal[NoDistance.NAN] | None = None + # the best known upper bound on distance + _distance_bound: int | Literal[NoDistance.NAN] | None = None _is_canonicalized: bool = False @@ -439,7 +442,7 @@ def get_code_params( bound: int | bool | None = None, decoder: decoders.DeferredErrorDecoderInput = None, **bound_kwargs: Any, - ) -> tuple[int, int, int | float]: + ) -> tuple[int, int, int | Literal[NoDistance.NAN]]: """Compute the parameters of this code: ``[n,k,d]``. Args: @@ -450,10 +453,11 @@ def get_code_params( **bound_kwargs: Keyword arguments to pass to get_distance_bound. Returns: - A tuple of integers, ``(n, k, d)``, where: + A tuple ``(n, k, d)``, where: - n is the number of data bits - k is the number of encoded ("logical") bits - - d is the code distance (or an upper bound on code distance). + - d is the code distance (or an upper bound on code distance), an integer, or + ``NO_DISTANCE`` if the code encodes nothing. """ dimension = self.dimension distance = self.get_distance(bound=bound, vector=None, decoder=decoder, **bound_kwargs) @@ -467,7 +471,7 @@ def get_distance( decoder: decoders.ErrorDecoderInput = None, method: DistanceMethod = "brouwer_zimmermann", **bound_kwargs: Any, - ) -> int | float: + ) -> int | Literal[NoDistance.NAN]: """Compute (or upper bound) the minimum Hamming weight of nontrivial code words. Args: @@ -487,7 +491,7 @@ def get_distance( **bound_kwargs: Keyword arguments to pass to get_distance_bound. Returns: - An integer distance (or bound) if it is defined, and np.nan otherwise. + An integer distance (or bound) if it is defined, and ``NO_DISTANCE`` otherwise. Raises: ValueError: If the method is unknown, or is ``"connected_support"`` for a nonbinary @@ -512,7 +516,7 @@ def get_distance_exact( vector: Sequence[int] | npt.NDArray[np.int_] | None = None, cutoff: int = 1, method: DistanceMethod = "brouwer_zimmermann", - ) -> int | float: + ) -> int | Literal[NoDistance.NAN]: """Compute the exact minimum Hamming weight of nontrivial code words. The cost of this computation generally grows exponentially with the size of the code. The @@ -535,7 +539,7 @@ def get_distance_exact( ``"connected_support"``. Returns: - An integer distance if it is defined, or np.nan otherwise. + An integer distance if it is defined, or ``NO_DISTANCE`` otherwise. Raises: ValueError: If the method is unknown, or is ``"connected_support"`` for a nonbinary @@ -547,6 +551,7 @@ def get_distance_exact( return known_distance # we do not know the exact distance, so compute it + distance: int | Literal[NoDistance.NAN] if vector is None and method == "connected_support": distance = _get_distance_connected_support( [(self.matrix, np.zeros((0, len(self)), dtype=np.uint8))], cutoff=cutoff @@ -609,7 +614,7 @@ def get_distance_exact( def get_distance_if_known( self, vector: Sequence[int] | npt.NDArray[np.int_] | None = None - ) -> int | float | None: + ) -> int | Literal[NoDistance.NAN] | None: """Retrieve a distance, or None if it is not known. Args: @@ -621,7 +626,7 @@ def get_distance_if_known( # the distance of dimension-0 codes is undefined if self.dimension == 0: - self._distance = np.nan + self._distance = NO_DISTANCE return self._distance @@ -634,7 +639,7 @@ def get_distance_bound( decoder: decoders.ErrorDecoderInput = None, backend: DistanceBackend = "auto", **decoder_kwargs: Any, # DEPRECATED: see qldpc._deprecations - ) -> int | float: + ) -> int | Literal[NoDistance.NAN]: """Use a randomized algorithm to compute an upper bound on code distance. The code distance is the minimum Hamming distance between two code words, or equivalently @@ -674,10 +679,10 @@ def get_distance_bound( ``decoder=`` instead. Returns: - An upper bound on distance if it is defined, or np.nan otherwise. If the exact distance - is already known, it is returned. When bounding code distance, the returned bound is - never worse than the best bound already known for this code, and the best bound is - remembered for later calls. + An upper bound on distance if it is defined, or ``NO_DISTANCE`` otherwise. If the exact + distance is already known, it is returned. When bounding code distance, the returned + bound is never worse than the best bound already known for this code, and the best bound + is remembered for later calls. Raises: TypeError: If num_trials is not an integer or None. @@ -1899,7 +1904,7 @@ def gauge_dimension(self) -> int: def get_code_params( self, *, bound: int | bool | None = None, **bound_kwargs: Any - ) -> tuple[int, int, int | float]: + ) -> tuple[int, int, int | Literal[NoDistance.NAN]]: """Compute the parameters of this code: ``[n,k,d]``. Args: @@ -1909,10 +1914,11 @@ def get_code_params( **bound_kwargs: Keyword arguments to pass to get_distance_bound. Returns: - A tuple of integers, ``(n, k, d)``, where: + A tuple ``(n, k, d)``, where: - n is the number of data qudits - k is the number of encoded ("logical") qudits - - d is the code distance (or an upper bound on code distance). + - d is the code distance (or an upper bound on code distance), an integer, or + ``NO_DISTANCE`` if the code encodes nothing. """ dimension = self.dimension distance = self.get_distance(bound=bound, **bound_kwargs) @@ -1924,7 +1930,7 @@ def get_distance( bound: int | bool | None = None, method: DistanceMethod = "brouwer_zimmermann", **bound_kwargs: Any, - ) -> int | float: + ) -> int | Literal[NoDistance.NAN]: """Compute (or upper bound) the minimum weight of nontrivial logical operators. Args: @@ -1940,7 +1946,7 @@ def get_distance( **bound_kwargs: Keyword arguments to pass to get_distance_bound. Returns: - An integer distance (or bound) if it is defined, and np.nan otherwise. + An integer distance (or bound) if it is defined, and ``NO_DISTANCE`` otherwise. Raises: ValueError: If the method is unknown, or is ``"connected_support"`` for a nonbinary @@ -1962,7 +1968,7 @@ def get_distance_exact( *, cutoff: int = 1, method: DistanceMethod = "brouwer_zimmermann", - ) -> int | float: + ) -> int | Literal[NoDistance.NAN]: """Compute the exact minimum weight of nontrivial logical operators. The weight of an operator is the number of qudits on which it acts nontrivially. For a @@ -1984,7 +1990,7 @@ def get_distance_exact( ``"connected_support"``. Returns: - An integer distance if it is defined, or np.nan otherwise. + An integer distance if it is defined, or ``NO_DISTANCE`` otherwise. Raises: ValueError: If the method is unknown, or is ``"connected_support"`` for a nonbinary @@ -1996,6 +2002,7 @@ def get_distance_exact( return known_distance # we do not know the exact distance, so compute it + distance: int | Literal[NoDistance.NAN] if method == "connected_support": # A dressed logical operator commutes with the stabilizers, and is trivial if it is in # the gauge group, which the rows of self.matrix generate. @@ -2042,14 +2049,14 @@ def get_distance_exact( self._distance = int(distance) return distance - def get_distance_if_known(self) -> int | float | None: + def get_distance_if_known(self) -> int | Literal[NoDistance.NAN] | None: """Retrieve a distance, if known. Otherwise, return None. """ # the distance of dimension-0 codes is undefined if self.dimension == 0: - self._distance = np.nan + self._distance = NO_DISTANCE return self._distance @@ -2061,7 +2068,7 @@ def get_distance_bound( cutoff: int | None = None, backend: DistanceBackend = "auto", **bound_kwargs: Any, - ) -> int | float: + ) -> int | Literal[NoDistance.NAN]: """Use a randomized algorithm to compute an upper bound on code distance. Specifically, use GAP's QDistRnd package to compute a distance bound. By default, @@ -2085,9 +2092,9 @@ def get_distance_bound( {gap_options}; see https://qec-pages.github.io/QDistRnd/doc/chap4.html. Returns: - An upper bound on distance if it is defined, or np.nan otherwise. The returned bound is - never worse than the best bound already known for this code, and the best bound is - remembered for later calls. + An upper bound on distance if it is defined, or ``NO_DISTANCE`` otherwise. The returned + bound is never worse than the best bound already known for this code, and the best bound + is remembered for later calls. Raises: NotImplementedError: If GAP cannot be called with ``backend="auto"``, or if the manual @@ -2758,10 +2765,10 @@ class CSSCode(QuditCode): _code_x: ClassicalCode _code_z: ClassicalCode - _distance_x: int | float | None = None - _distance_z: int | float | None = None - _distance_bound_x: int | float | None = None - _distance_bound_z: int | float | None = None + _distance_x: int | Literal[NoDistance.NAN] | None = None + _distance_z: int | Literal[NoDistance.NAN] | None = None + _distance_bound_x: int | Literal[NoDistance.NAN] | None = None + _distance_bound_z: int | Literal[NoDistance.NAN] | None = None _equal_distance_xz: bool # are the X and Z distances promised to be equal? @@ -3480,7 +3487,7 @@ def get_code_params( bound: int | bool | None = None, decoder: decoders.DeferredErrorDecoderInput = None, **bound_kwargs: Any, - ) -> tuple[int, int, int | float]: + ) -> tuple[int, int, int | Literal[NoDistance.NAN]]: """Compute the parameters of this code: ``[n,k,d]``. Args: @@ -3491,10 +3498,11 @@ def get_code_params( **bound_kwargs: Keyword arguments to pass to get_distance_bound. Returns: - A tuple of integers, ``(n, k, d)``, where: + A tuple ``(n, k, d)``, where: - n is the number of data qudits - k is the number of encoded ("logical") qudits - - d is the code distance (or an upper bound on code distance). + - d is the code distance (or an upper bound on code distance), an integer, or + ``NO_DISTANCE`` if the code encodes nothing. """ dimension = self.dimension distance = self.get_distance(bound=bound, decoder=decoder, **bound_kwargs) @@ -3508,7 +3516,7 @@ def get_distance( decoder: decoders.DeferredErrorDecoderInput = None, method: DistanceMethod = "brouwer_zimmermann", **bound_kwargs: Any, - ) -> int | float: + ) -> int | Literal[NoDistance.NAN]: """Compute (or upper bound) the minimum weight of nontrivial logical operators. Args: @@ -3529,7 +3537,7 @@ def get_distance( **bound_kwargs: Keyword arguments to pass to get_distance_bound. Returns: - An integer distance (or bound) if it is defined, and np.nan otherwise. + An integer distance (or bound) if it is defined, and ``NO_DISTANCE`` otherwise. Raises: ValueError: If the method is unknown, or is ``"connected_support"`` for a nonbinary @@ -3554,7 +3562,7 @@ def get_distance_exact( *, cutoff: int = 1, method: DistanceMethod = "brouwer_zimmermann", - ) -> int | float: + ) -> int | Literal[NoDistance.NAN]: """Compute the exact minimum weight of nontrivial logical operators. For a subsystem code, logical operators are dressed: they may be multiplied by gauge @@ -3580,7 +3588,7 @@ def get_distance_exact( enumeration, and reject ``"connected_support"``. Returns: - An integer distance if it is defined, or np.nan otherwise. + An integer distance if it is defined, or ``NO_DISTANCE`` otherwise. Raises: ValueError: If the method is unknown, or is ``"connected_support"`` for a nonbinary @@ -3677,7 +3685,9 @@ def get_distance_exact( return distance - def get_distance_if_known(self, pauli: PauliXZLike | None = None) -> int | float | None: + def get_distance_if_known( + self, pauli: PauliXZLike | None = None + ) -> int | Literal[NoDistance.NAN] | None: """Retrieve a distance, if known. Otherwise, return None. @@ -3686,7 +3696,7 @@ def get_distance_if_known(self, pauli: PauliXZLike | None = None) -> int | float # the distances of dimension-0 codes are undefined if self.dimension == 0: - self._distance = self._distance_x = self._distance_z = np.nan + self._distance = self._distance_x = self._distance_z = NO_DISTANCE if pauli is Pauli.X: return self._distance_x @@ -3711,7 +3721,7 @@ def get_distance_bound( decoder: decoders.DeferredErrorDecoderInput = None, backend: DistanceBackend = "auto", **bound_kwargs: Any, - ) -> int | float: + ) -> int | Literal[NoDistance.NAN]: """Estimate an upper bound on code distance with a randomized algorithm. Set ``backend="gap"`` to use GAP's QDistRnd package or ``backend="sqetch"`` to use the @@ -3745,9 +3755,10 @@ def get_distance_bound( that backend, and raise a ValueError if it is unavailable. Returns: - An upper bound on distance if it is defined, or np.nan otherwise. If the exact distance - is already known, it is returned. The returned bound is never worse than the best bound - already known for this code, and the best bound is remembered for later calls. + An upper bound on distance if it is defined, or ``NO_DISTANCE`` otherwise. If the exact + distance is already known, it is returned. The returned bound is never worse than the + best bound already known for this code, and the best bound is remembered for later + calls. Raises: NotImplementedError: If the manual GAP workflow is declined with ``backend="gap"``. @@ -3843,7 +3854,7 @@ def get_distance_bound_with_decoder( cutoff: int | None = None, decoder: decoders.DeferredErrorDecoderInput = None, **decoder_kwargs: Any, # DEPRECATED: see qldpc._deprecations - ) -> int | float: + ) -> int | Literal[NoDistance.NAN]: r"""Use a randomized algorithm to compute an upper bound on code distance. Specifically, use the algorithm described in arXiv:2308.07915, also explained below. @@ -3899,6 +3910,8 @@ def get_distance_bound_with_decoder( define equivalent decoding problems, so each problem is decoded at most once. At most 100 problems beyond ``num_trials`` are drawn to replace problems on which the decoder fails. + A code without logical qudits has no distance, and gives ``NO_DISTANCE``. + Raises: ValueError: If the decoder fails to infer an error that reproduces its syndrome in every decoding problem that it attempts. @@ -3906,6 +3919,8 @@ def get_distance_bound_with_decoder( pauli = Pauli.coerce_xz(pauli) reject_prebuilt_decoder(decoder, _CSS_DISTANCE_BOUND_REASON) cutoff = cutoff or 0 + if self.dimension == 0: + return NO_DISTANCE # pretend without loss of generality that we are computing the X-distance pauli_z: PauliXZ = Pauli.Z if pauli is Pauli.X else Pauli.X @@ -4369,7 +4384,9 @@ def _get_distance_of_sectors( get_automorphisms=lambda: self._found_qudit_automorphisms, ) - def _get_distance_exact(self, pauli: PauliXZ | None, *, method: DistanceMethod) -> int | float: + def _get_distance_exact( + self, pauli: PauliXZ | None, *, method: DistanceMethod + ) -> int | Literal[NoDistance.NAN]: """Method for subclasses to compute specialized exact distance calculations. Args: @@ -4378,7 +4395,9 @@ def _get_distance_exact(self, pauli: PauliXZ | None, *, method: DistanceMethod) """ return NotImplemented - def _get_distance_bound_if_known(self, pauli: PauliXZ | None = None) -> int | float | None: + def _get_distance_bound_if_known( + self, pauli: PauliXZ | None = None + ) -> int | Literal[NoDistance.NAN] | None: """Retrieve the best known upper bound on (X-, Z-, or ordinary) distance, or ``None``.""" if pauli is Pauli.X: candidates = [self._distance_x, self._distance_bound_x] @@ -4395,7 +4414,9 @@ def _get_distance_bound_if_known(self, pauli: PauliXZ | None = None) -> int | fl ] return min((bound for bound in candidates if bound is not None), default=None) - def _update_distance_bound(self, bound: float, pauli: PauliXZ) -> float: + def _update_distance_bound( + self, bound: int | Literal[NoDistance.NAN], pauli: PauliXZ + ) -> int | Literal[NoDistance.NAN]: """Remember an upper bound on the X-distance or Z-distance, and return the best bound.""" known_bound = self._get_distance_bound_if_known(None if self._equal_distance_xz else pauli) bound = bound if known_bound is None else min(bound, known_bound) diff --git a/src/qldpc/codes/common_test.py b/src/qldpc/codes/common_test.py index 94c50afdb..85332aef5 100644 --- a/src/qldpc/codes/common_test.py +++ b/src/qldpc/codes/common_test.py @@ -211,8 +211,8 @@ def test_distance_classical(bits: int = 3) -> None: trivial_code = codes.ClassicalCode([[1, 0], [1, 1]]) random_vector = np.random.randint(2, size=len(trivial_code)) assert trivial_code.dimension == 0 - assert np.isnan(trivial_code.get_distance_exact()) - assert np.isnan(trivial_code.get_distance_bound()) + assert trivial_code.get_distance_exact() is codes.NO_DISTANCE + assert trivial_code.get_distance_bound() is codes.NO_DISTANCE assert ( np.count_nonzero(random_vector) == trivial_code.get_distance_exact(vector=random_vector) @@ -591,7 +591,7 @@ def test_distance_qudit() -> None: assert code.get_distance(bound=True) == -1 # the distance of dimension-0 codes is undefined - assert np.isnan(codes.QuditCode([[0, 1]]).get_distance()) + assert codes.QuditCode([[0, 1]]).get_distance() is codes.NO_DISTANCE # fallback pythonic brute-force distance calculation code = codes.QuditCode(codes.SurfaceCode(2, field=3).matrix) @@ -1711,8 +1711,10 @@ def test_distance_css() -> None: trivial_code = codes.ClassicalCode([[1, 0], [1, 1]]) code = codes.HGPCode(trivial_code) assert code.dimension == 0 - assert np.isnan(code.get_distance(bound=True)) - assert np.isnan(code.get_distance(bound=False)) + assert code.get_distance(bound=True) is codes.NO_DISTANCE + assert code.get_distance(bound=False) is codes.NO_DISTANCE + with unittest.mock.patch("qldpc.codes.common.build_error_decoder", side_effect=AssertionError): + assert code.get_distance_bound_with_decoder(Pauli.X) is codes.NO_DISTANCE def test_css_decoder_distance_bound_skips_gap_probe() -> None: diff --git a/src/qldpc/codes/distance/README.md b/src/qldpc/codes/distance/README.md index 1cbe89d7a..b8be4c5e2 100644 --- a/src/qldpc/codes/distance/README.md +++ b/src/qldpc/codes/distance/README.md @@ -10,7 +10,7 @@ Users reach it through the `get_distance_exact` method of binary codes in `src/q `__init__.py` re-exports these names, and everything else in this directory is private: -- `get_distance_classical`, `get_distance_quantum`, `DistanceMethod`, and `DistanceBackend`, from `core.py`. +- `get_distance_classical`, `get_distance_quantum`, `DistanceMethod`, `DistanceBackend`, and `NoDistance`, whose only member `NoDistance.NAN` is also exported as `NO_DISTANCE`, from `core.py`. - `CSSLogicalCounter`, `CSSLogicalCounts`, `count_logical_ops_css`, `LogicalCountMethod`, and `compile_search_setup`, from `logical_counts.py`. ## Terms @@ -21,7 +21,8 @@ Users reach it through the `get_distance_exact` method of binary codes in `src/q - **Information set:** in the code, the basis of a sector, row reduced so that a block of its columns forms an identity matrix. `brouwer_zimmermann.md` calls this basis a search basis, and uses "information set" for the block of columns itself. - **Basis family** (`BasisFamily`): the information sets that one Brouwer–Zimmermann search or count uses together, stacked into arrays. -- **Certificate** (`Certificate`): a lower bound on the weight of the words that a Brouwer–Zimmermann search has not yet visited, which tells the search when it can stop. +- **Certificate** (`Certificate`): a lower bound on the weight of the words that a Brouwer–Zimmermann search has neither visited nor skipped, which tells the search when it can stop. + With completion tables, it bounds only the words lighter than the best weight (see Section 3.9 of `brouwer_zimmermann.md`). - **Connected-support search** (`SupportSearch`): the depth-first search of the supports that grow through unsatisfied parity checks, for a binary matrix of checks `H` and an excluded matrix `G`, whose words are binary, or symplectic for qubit codes that are not CSS. It is another exact method, selected by `method="connected_support"`, and it also finds light vectors with a given syndrome. - **Scan:** the enumeration of the combinations of rows of each information set, layer by layer. diff --git a/src/qldpc/codes/distance/__init__.py b/src/qldpc/codes/distance/__init__.py index 8a57867ea..051545a53 100644 --- a/src/qldpc/codes/distance/__init__.py +++ b/src/qldpc/codes/distance/__init__.py @@ -33,8 +33,10 @@ from qldpc import _deprecations from .core import ( + NO_DISTANCE, DistanceBackend, DistanceMethod, + NoDistance, get_distance_classical, get_distance_quantum, ) @@ -47,11 +49,13 @@ ) __all__ = [ + "NO_DISTANCE", "CSSLogicalCounter", "CSSLogicalCounts", "DistanceBackend", "DistanceMethod", "LogicalCountMethod", + "NoDistance", "compile_search_setup", "count_logical_ops_css", "get_distance_classical", diff --git a/src/qldpc/codes/distance/_binary.py b/src/qldpc/codes/distance/_binary.py index 5af4b6ef0..95b16640b 100644 --- a/src/qldpc/codes/distance/_binary.py +++ b/src/qldpc/codes/distance/_binary.py @@ -62,22 +62,28 @@ def assert_binary(vectors: npt.ArrayLike, name: str) -> None: def as_binary_matrix( - vectors: npt.ArrayLike, *, num_cols: int | None = None + vectors: npt.ArrayLike, name: str = "rows", *, num_cols: int | None = None ) -> npt.NDArray[np.uint8]: - """Convert a collection of binary rows to a two-dimensional uint8 array.""" + """Convert a collection of binary rows to a two-dimensional uint8 array. + + The dimension and width checks use ``name`` in their error messages, as ``assert_binary`` does. + If NumPy cannot convert the rows to an array, such as when their lengths differ, it raises its + own error. + """ matrix = np.asarray(vectors, dtype=np.uint8) if matrix.size == 0: if num_cols is None: if matrix.ndim != 2: - raise ValueError("Cannot infer the width of an empty binary matrix") + raise ValueError(f"Cannot infer the width of {name}, which are empty") num_cols = matrix.shape[1] return np.empty((0, num_cols), dtype=np.uint8) matrix = np.atleast_2d(matrix) if matrix.ndim != 2: - raise ValueError("Binary generators must be a two-dimensional matrix") + raise ValueError(f"The {name} must be a two-dimensional matrix") if num_cols is not None and matrix.shape[1] != num_cols: raise ValueError( - f"Binary matrices have incompatible widths {num_cols} and {matrix.shape[1]}" + f"The {name} have {matrix.shape[1]} columns, which is incompatible with the width" + f" {num_cols}" ) return matrix @@ -86,7 +92,15 @@ def get_nested_code_basis( logical_ops: npt.NDArray[np.uint8], stabilizers: npt.NDArray[np.uint8], ) -> tuple[npt.NDArray[np.uint8], npt.NDArray[np.uint8]]: - """Build an adapted basis for ``span(stabilizers, logical_ops) / span(stabilizers)``.""" + """Build an adapted basis for ``span(stabilizers, logical_ops) / span(stabilizers)``. + + Returns: + An independent basis whose first rows span the stabilizers, and the quotient labels of its + rows: row ``i`` of the labels holds the coordinates of row ``i`` of the basis in the + quotient, which are zero exactly for the stabilizer rows. If the quotient is trivial, the + labels have no columns. Callers must handle that case first, since ``get_eligible_rows`` + and the scans read labels with no columns as excluding nothing. + """ stabilizer_basis = get_independent_rows(stabilizers) stabilizer_pivots = [int(np.flatnonzero(row)[0]) for row in stabilizer_basis] reduced_logical_ops = logical_ops.copy() diff --git a/src/qldpc/codes/distance/_binary_test.py b/src/qldpc/codes/distance/_binary_test.py index dcc40dbb2..b15e99c91 100644 --- a/src/qldpc/codes/distance/_binary_test.py +++ b/src/qldpc/codes/distance/_binary_test.py @@ -25,10 +25,10 @@ def test_binary_matrices_and_nested_bases() -> None: assert binary.as_binary_matrix([], num_cols=3).shape == (0, 3) with pytest.raises(ValueError, match="infer the width"): binary.as_binary_matrix([]) - with pytest.raises(ValueError, match="two-dimensional"): + with pytest.raises(ValueError, match="The rows must be a two-dimensional"): binary.as_binary_matrix(np.zeros((1, 1, 1))) - with pytest.raises(ValueError, match="incompatible widths"): - binary.as_binary_matrix([[1, 0]], num_cols=3) + with pytest.raises(ValueError, match="The stabilizers have 2 columns"): + binary.as_binary_matrix([[1, 0]], "stabilizers", num_cols=3) matrix = np.array([[0, 1, 1], [1, 1, 0], [1, 0, 1], [0, 0, 0]], dtype=np.uint8) basis, tags = binary.get_nested_code_basis(matrix, np.array([[1, 0, 1]], dtype=np.uint8)) diff --git a/src/qldpc/codes/distance/_counting/scan.py b/src/qldpc/codes/distance/_counting/scan.py index 86522668d..25aea8669 100644 --- a/src/qldpc/codes/distance/_counting/scan.py +++ b/src/qldpc/codes/distance/_counting/scan.py @@ -81,7 +81,7 @@ def count_words( labels, weight, get_disjoint_basis_family(), - *information_sets.get_overlapping_family_policy(tuning, len(basis)), + information_sets.get_overlapping_family_policy(tuning, len(basis)), ) # Matroid partitioning runs outside of _select_basis_family, which can be compiled. partitioned = information_sets.get_partitioned_family( @@ -190,21 +190,17 @@ def _select_basis_family( labels: npt.NDArray[np.uint8], weight: int, basis_family: information_sets.BasisFamily, - max_sets: int, - minimum_work: int, - attempts_per_set: int, - binomials: npt.NDArray[np.int64], - rng: np.random.Generator, + policy: information_sets.OverlappingFamilyPolicy, ) -> information_sets.BasisFamily: """Run the basis family stages that can be compiled: pruning, then an overlapping replacement. The stages aim at the counting stop ``L(t) >= weight + 1``; any basis family gives an exact - count. + count. ``policy`` is that of the overlapping replacement. """ upper_bound = weight + 1 pruned = information_sets.prune_family(basis_family, upper_bound) replacement = information_sets.get_overlapping_family( - basis, labels, pruned, upper_bound, max_sets, minimum_work, attempts_per_set, binomials, rng + basis, labels, pruned, upper_bound, policy ) return replacement if len(replacement.ranks) > 0 else pruned diff --git a/src/qldpc/codes/distance/_search/affine_search.py b/src/qldpc/codes/distance/_search/affine_search.py index 1def129df..ba6f968a7 100644 --- a/src/qldpc/codes/distance/_search/affine_search.py +++ b/src/qldpc/codes/distance/_search/affine_search.py @@ -101,7 +101,6 @@ def scan_family( record=False, css_weight=False, chunked=True, - threshold=False, ranked=False, record_support=False, ) @@ -143,8 +142,9 @@ def scan_information_set( numbers of supports that the chunk weighed, of bounds that it looked up, and of supports that it skipped. """ - # This scan neither reads ranks nor weighs X parts nor records words or supports, nor takes a - # skip bound. + # This scan neither reads ranks nor weighs X parts nor records words or supports. Its skip + # bound is the best weight at the start, which the best weight never exceeds, so it skips + # exactly the completions that cannot lower the best weight. return _scan( generators, labels, @@ -162,7 +162,7 @@ def scan_information_set( block_bits, prefixes[0, 0], prefixes[0, 0], - 0, + best, support, one_set=True, pair_steps=True, @@ -171,7 +171,6 @@ def scan_information_set( record=False, css_weight=False, chunked=True, - threshold=False, ranked=False, record_support=False, ) @@ -236,7 +235,6 @@ def scan_coset_layer( record=True, css_weight=False, chunked=False, - threshold=False, ranked=False, record_support=False, ) @@ -290,7 +288,6 @@ def scan_coset_layer_by_css_weight( record=True, css_weight=True, chunked=False, - threshold=False, ranked=False, record_support=False, ) @@ -355,7 +352,6 @@ def scan_coset_layer_by_tails( record=False, css_weight=False, chunked=False, - threshold=False, ranked=True, record_support=True, ) @@ -392,16 +388,22 @@ def scan_coset_layer_by_tails_with_tables( ``rank``, and none for the rows of the prefix at or above ``rank``. Its last two rows have at least as many pivot ones as the last two rows of the basis, the pair with the fewest pivot rows, since the pivot rows come first: ``max(0, 2 - (num_rows - rank))``. Its tail weighs at least - the bound of the tables. If the sum reaches ``skip_bound``, the scan skips these words. + the bound of the tables. If the sum reaches the lesser of ``skip_bound`` and the best weight so + far, which starts at ``best``, the scan skips these words. A coset query passes as ``skip_bound`` the least bound that proves a word of its coset heavier than its maximum weight: one more than the greatest Hamming weight, at most ``cutoff``, that a - word of the coset can have. So the scan skips no word within that maximum weight. - ``skip_bound`` must be at most ``cutoff + 1``. Then, while the best weight exceeds ``cutoff``, - the scan skips every completion whose bound shows that it cannot lower the best weight, and - more. A larger skip bound would make it weigh completions that cannot lower the best weight. - Skipped words can be lighter than the best weight, which is then the weight of the lightest word - weighed. + word of the coset can have, which is at most ``cutoff + 1``. The scan returns once the best + weight is at most ``cutoff``, so the best weight reaches the skip bound at every lookup, and the + scan skips at the skip bound. So it skips no word within that maximum weight, and every + completion whose bound shows that it cannot lower the best weight, and more. Skipped words can + be lighter than the best weight, which is then the weight of the lightest word weighed. + + A search for the minimum weight of a coset passes a skip bound of at least ``best``, so the scan + skips, as a distance search does, the completions whose bound reaches the best weight so far. + So it skips no word lighter than the best weight at the time, and finds the lightest word of the + layer if it is lighter than ``best``, unless it returns at the cutoff first (see + ``coset_leader_search._search_coset_for_min_weight``). Returns: The weight of the lightest word found, or ``best`` if none is lighter, and the numbers of @@ -436,74 +438,6 @@ def scan_coset_layer_by_tails_with_tables( record=False, css_weight=False, chunked=False, - threshold=True, - ranked=True, - record_support=True, - ) - return best, weighed, skipped - - -@always_compiled -def scan_coset_layer_by_tails_with_tables_below_best( - offset_tail: npt.NDArray[np.uint64], - tails: npt.NDArray[np.uint64], - ranks: npt.NDArray[np.int64], - widths: npt.NDArray[np.int64], - support: npt.NDArray[np.int64], - prefixes: npt.NDArray[np.uint64], - best: int, - cutoff: int, - tables: npt.NDArray[np.uint8], - offsets: npt.NDArray[np.int64], - block_bits: int, - best_support: npt.NDArray[np.int64], -) -> tuple[int, int, int]: - """Scan a layer by tails, and skip by tables the completions that cannot lower the best weight. - - The scan weighs the words of the layer as ``scan_coset_layer_by_tails_with_tables`` does, but - skips, as a distance search does, the completions whose bound reaches the best weight so far, - which starts at ``best``, rather than those whose bound reaches a fixed skip bound. So it skips - no word lighter than the best weight at the time, and finds the lightest word of the layer if it - is lighter than ``best``, unless it returns at the cutoff first. A search for the minimum - weight of a coset passes as ``best`` the skip bound of its maximum weight until it finds a word - within that weight, and then the Hamming weight of its lightest word, or one less where every - word of the coset has its parity (see ``coset_leader_search._search_coset_for_min_weight``). - - Returns: - The weight of the lightest word found, or ``best`` if none is lighter, and the numbers of - words weighed and of supports skipped. - """ - prefixes[0, 0] = offset_tail - support[0] = 0 - best, _, weighed, _, skipped = _scan( - tails, - tails, - ranks, - widths, - widths, - support, - prefixes, - 0, - best, - cutoff, - 0, - tables, - offsets, - block_bits, - # This scan neither reads labels nor weighs X parts nor records words, nor takes a skip - # bound. - prefixes[0, 0], - prefixes[0, 0], - 0, - best_support, - one_set=True, - pair_steps=True, - deleted_pivots=True, - test_labels=False, - record=False, - css_weight=False, - chunked=False, - threshold=False, ranked=True, record_support=True, ) @@ -541,30 +475,27 @@ def _scan( record: bool, css_weight: bool, chunked: bool, - threshold: bool, ranked: bool, record_support: bool, ) -> tuple[int, int, int, int, int]: """Scan the supports of a layer of a stack of search bases, starting from ``support``. The arguments are those of ``scan_family``, ``scan_information_set``, - ``scan_coset_layer_by_css_weight``, ``scan_coset_layer_by_tails``, - ``scan_coset_layer_by_tails_with_tables``, and - ``scan_coset_layer_by_tails_with_tables_below_best``, and the flags, which callers pass as - constants, choose what the scan does: + ``scan_coset_layer_by_css_weight``, ``scan_coset_layer_by_tails``, and + ``scan_coset_layer_by_tails_with_tables``, and the flags, which callers pass as constants, + choose what the scan does: - ``one_set``: scan only the first basis, rather than all bases in lockstep. - ``pair_steps``: with ``one_set``, step through the first of the last two rows of a support - separately, and skip completions by the completion tables, if any, as - ``scan_information_set`` describes. + separately, and skip the completions whose bound by the completion tables, if any, reaches the + lesser of ``skip_bound`` and the best weight so far (see ``scan_information_set`` and + ``scan_coset_layer_by_tails_with_tables``). - ``deleted_pivots``: the rows have their pivot columns deleted, so a word also weighs the number of its rows whose pivots are deleted, which is ``len(support)`` for an information set. - ``test_labels``: keep only words whose quotient labels make them eligible. - ``record``: copy each lighter word into ``representative``. - ``css_weight``: weigh words by their CSS weight, with ``x_mask``, rather than Hamming weight. - ``chunked``: return once the weighed supports and the lookups reach ``batch_size`` together. - - ``threshold``: with ``pair_steps``, skip completions whose bound reaches ``skip_bound``, - rather than ``best`` (see ``scan_coset_layer_by_tails_with_tables``). - ``ranked``: with ``one_set`` and ``deleted_pivots``, lift the assumption of a scan of one basis that the basis is an information set: only its first ``ranks[0]`` rows hold pivots. A word then weighs the number of its rows below ``ranks[0]``, rather than ``len(support)``, plus @@ -599,10 +530,10 @@ def _scan( return best, depth, evaluated, lookups, skipped first = support[pair_depth] - # Skip every completion from row first on if none can be lighter than the best word, or - # with a threshold, if every one weighs at least the skip bound. The check for tables - # runs once for each prefix: as a constant, it would let numba merge the lookup into the - # loop, and then count references to the tables for each lookup. + # Skip every completion from row first on if every one weighs at least the lesser of the + # skip bound and the best weight. The check for tables runs once for each prefix: as a + # constant, it would let numba merge the lookup into the loop, and then count references + # to the tables for each lookup. prefix = prefixes[pair_depth, 0] if len(tables): lookups += 1 @@ -617,7 +548,7 @@ def _scan( for selected in range(pair_depth): if support[selected] < ranks[0]: pivot_ones += 1 - limit = (skip_bound if threshold else best) - pivot_ones + limit = min(skip_bound, best) - pivot_ones if get_completion_bound(tables, offsets, block_bits, first, prefix, limit) >= limit: remaining = num_rows - first skipped += remaining * (remaining - 1) // 2 diff --git a/src/qldpc/codes/distance/_search/affine_search_test.py b/src/qldpc/codes/distance/_search/affine_search_test.py index df636a791..ed48628ab 100644 --- a/src/qldpc/codes/distance/_search/affine_search_test.py +++ b/src/qldpc/codes/distance/_search/affine_search_test.py @@ -473,11 +473,13 @@ def test_coset_tail_scans_below_the_best_skip_only_completions_that_cannot_lower ) -> None: """A coset tail scan below the best weight skips completions only if they cannot lower it. - The bases of six rows have rank deficiency zero and three. From a best weight of at most the - least weight of a layer, which no word lowers, the scan skips the supports that a scan with that - skip bound skips, which the test above checks against enumeration. From a best weight above - every word of the layer, it finds a lightest word, since every word that it skips weighs at - least the best weight at the time. It weighs or skips every support, and skips some. + The scan takes its best weight at the start as its skip bound, as a search for the minimum + weight of a coset does. The bases of six rows have rank deficiency zero and three. From a best + weight of at most the least weight of a layer, which no word lowers, the scan skips the supports + that a scan with that skip bound skips, which the test above checks against enumeration. From a + best weight above every word of the layer, it finds a lightest word, since every word that it + skips weighs at least the best weight at the time. It weighs or skips every support, and skips + some. """ rng = np.random.default_rng(pytestconfig.getoption("randomly_seed")) num_rows, num_columns, block_bits = 6, 24, 4 @@ -504,10 +506,8 @@ def test_coset_tail_scans_below_the_best_skip_only_completions_that_cannot_lower np.empty((layer, 1, widths[0]), dtype=np.uint64), ) best_support = np.full(layer, -1, dtype=np.int64) - found, weighed, skipped = ( - affine_search.scan_coset_layer_by_tails_with_tables_below_best( - *arguments, best, -1, *table_arguments, best_support - ) + found, weighed, skipped = affine_search.scan_coset_layer_by_tails_with_tables( + *arguments, best, -1, best, *table_arguments, best_support ) assert weighed + skipped == math.comb(num_rows, layer) if best > least: diff --git a/src/qldpc/codes/distance/_search/brouwer_zimmermann.py b/src/qldpc/codes/distance/_search/brouwer_zimmermann.py index 5e415043a..07c62677e 100644 --- a/src/qldpc/codes/distance/_search/brouwer_zimmermann.py +++ b/src/qldpc/codes/distance/_search/brouwer_zimmermann.py @@ -56,27 +56,50 @@ def get_distance_brouwer_zimmermann_many( get_automorphisms: Callable[[], npt.NDArray[np.int_] | None] | None = None, tuning: SearchTuning | None = None, ) -> int: - """Advance nested-code searches round-robin while sharing their best upper bound. + """Compute the minimum distance across nested-code problems by BZ searches. - ``get_disjoint_basis_families``, if given, holds a function for each problem that returns a - disjoint basis family of it: sets with pairwise disjoint pivot blocks that span the problem's - basis, with quotient labels of the problem, certified with this ``weight_divisor``. The result - of ``information_sets.get_disjoint_family`` is one, and a caller that keeps such families saves - that first stage. Each function is called only if its search needs it. + The searches advance round-robin, and share their best upper bound: the weight of the lightest + witness that any of them has found. - ``get_automorphisms``, if given, returns an array of generators of a group of column - permutations, or None. Row ``i`` of the array maps column ``j`` to column ``row[j]``. Each - generator must map the searched space and the excluded subspace of every problem onto - themselves. A search calls the function before a costly layer, verifies the generators for its - problem, and may then scan fewer information sets (see ``automorphism_families``). + Args: + problems: Nested-code bases and their quotient labels (see ``get_nested_code_basis``). The + distance of a problem is the least weight of a word of its basis with nonzero labels, or + of any nonzero word if its labels have no columns. Each problem must have such a word. + cutoff: Return once the best upper bound falls to or below this cutoff. + weight_divisor: A proven divisor of the weight of every word whose distance a problem + minimizes over, with which the certificates round their lower bounds up. Default: 1. + upper_bound: A cap on the result, if any, such as the weight of a known witness that is + not in ``problems``, or ``max_weight + 1`` under a weight limit. The searches stop + once they certify that every word is at least this heavy. + get_disjoint_basis_families: If given, holds a function for each problem that returns a + disjoint basis family of it: sets with pairwise disjoint pivot blocks that span the + problem's basis, with quotient labels of the problem, certified with this + ``weight_divisor``. The result of ``information_sets.get_disjoint_family`` is one, + and a caller that keeps such families saves that first stage. Each function is called + only if its search needs it. + get_automorphisms: If given, returns an array of generators of a group of column + permutations, or None. Row ``i`` of the array maps column ``j`` to column ``row[j]``. + Each generator must map the searched space and the excluded subspace of every problem + onto themselves. A search calls the function before a costly layer, verifies the + generators for its problem, and may then scan fewer information sets (see + ``automorphism_families``). + tuning: Sets the thresholds that choose the route of each search (see ``_tuning``), and + defaults to the active tuning, which every search then uses throughout. - ``tuning`` sets the thresholds that choose the route of each search (see ``_tuning``), and - defaults to the active tuning, which every search then uses throughout. + Returns: + The minimum of the distances of the problems and ``upper_bound``, or, after an early exit + at ``cutoff``, the weight of a witness or ``upper_bound``, whichever is at most + ``cutoff``. Raises: - ValueError: If a generator from ``get_automorphisms`` is not an automorphism of a problem - whose search uses it. + ValueError: If there are neither problems nor an upper bound, or if a generator from + ``get_automorphisms`` is not an automorphism of a problem whose search uses it. """ + if not problems: + if upper_bound is None: + raise ValueError("Require at least one problem or an upper bound") + return upper_bound + tuning = get_tuning() if tuning is None else tuning initial_upper_bounds = [ _get_lightest_eligible_row_weight(basis, labels) for basis, labels in problems @@ -131,7 +154,7 @@ def get_distance_brouwer_zimmermann( upper_bound: int | None = None, tuning: SearchTuning | None = None, ) -> int: - """Compute an exact nested-code distance with the Brouwer-Zimmermann algorithm. + """Compute a nested-code distance with the Brouwer-Zimmermann algorithm. The disjoint and rank-deficient bound follows Algorithm 1 of https://arxiv.org/abs/1603.06757. The overlapping-set extension divides pivot occupancy by the maximum coordinate coverage. @@ -141,7 +164,9 @@ def get_distance_brouwer_zimmermann( An ``upper_bound`` starts the search as though a witness of that weight were known, so the search returns ``min(distance, upper_bound)`` and stops once it certifies that no word is - lighter than ``upper_bound``. ``tuning`` defaults to the active tuning. + lighter than ``upper_bound``. ``tuning`` defaults to the active tuning. As with + ``get_distance_brouwer_zimmermann_many``, the basis needs a nontrivial word, and an early exit + once a witness weighs at most ``cutoff`` returns its weight, which can exceed the distance. """ return get_distance_brouwer_zimmermann_many( [(basis, labels)], @@ -595,7 +620,7 @@ def _prepare_brouwer_zimmermann_search( labels, family, best, - *information_sets.get_overlapping_family_policy(tuning, dimension), + information_sets.get_overlapping_family_policy(tuning, dimension), ), lambda family, best: information_sets.get_partitioned_family( basis, diff --git a/src/qldpc/codes/distance/_search/brouwer_zimmermann_test.py b/src/qldpc/codes/distance/_search/brouwer_zimmermann_test.py index 5ad42a483..0d356c5e1 100644 --- a/src/qldpc/codes/distance/_search/brouwer_zimmermann_test.py +++ b/src/qldpc/codes/distance/_search/brouwer_zimmermann_test.py @@ -844,6 +844,7 @@ def assert_skips(distance: int) -> int: homogeneous=True, method="brute_force", ) + assert expected is not qldpc.codes.NO_DISTANCE for upper_bound in (None, expected + 1): skipped = completion_tables.skipped # One sixteen-column block covers the tails of this code, which makes its bounds exact, @@ -879,10 +880,11 @@ def assert_skips(distance: int) -> int: five_qubit_code = qldpc.codes.FiveQubitCode() logicals, stabilizers = five_qubit_code.get_logical_ops(), five_qubit_code.get_stabilizer_ops() skipped = completion_tables.skipped - distance = qldpc.codes.get_distance_quantum( + quantum_distance = qldpc.codes.get_distance_quantum( logicals, stabilizers, homogeneous=False, cutoff=0, method="brouwer_zimmermann" ) - assert assert_skips(distance) == qldpc.codes.get_distance_quantum( + assert quantum_distance is not qldpc.codes.NO_DISTANCE + assert assert_skips(quantum_distance) == qldpc.codes.get_distance_quantum( logicals, stabilizers, homogeneous=False, method="brute_force" ) @@ -1196,3 +1198,10 @@ def recording_select(*args: Any, **kwargs: Any) -> completion_bounds.CompletionT get_automorphisms=code._get_qudit_automorphisms, ) == code.get_distance(method="brute_force") assert block_bits == expected_block_bits + + +def test_many_problems_need_a_problem_or_an_upper_bound() -> None: + """Without problems, the result is the known upper bound, which must then be given.""" + with pytest.raises(ValueError, match="at least one problem or an upper bound"): + bz.get_distance_brouwer_zimmermann_many([], cutoff=0) + assert bz.get_distance_brouwer_zimmermann_many([], cutoff=0, upper_bound=3) == 3 diff --git a/src/qldpc/codes/distance/_search/brute_force.py b/src/qldpc/codes/distance/_search/brute_force.py index 45a8ae791..f967f98ed 100644 --- a/src/qldpc/codes/distance/_search/brute_force.py +++ b/src/qldpc/codes/distance/_search/brute_force.py @@ -27,7 +27,7 @@ def get_distance_quantum_brute_force( Every combination of the rows of logical_ops that uses at least one row, multiplied by every combination of the rows of stabilizers, is a candidate. The rows of logical_ops and stabilizers - must together be linearly independent; callers in this package construct or verify such bases. + must together be linearly independent; callers in this package construct such bases. Args: logical_ops: Binary rows that, together with the stabilizers, generate the logical space. @@ -40,7 +40,9 @@ def get_distance_quantum_brute_force( Returns: The minimum weight of a nontrivial logical operator or, after an early exit, the weight (at - most cutoff) of the lightest operator seen. + most cutoff) of the lightest operator seen. Without logical_ops, the result is one more + than the largest possible weight: the number of columns if homogeneous, and otherwise the + number of qubits. """ num_bits = np.shape(logical_ops)[-1] @@ -85,11 +87,17 @@ def get_distance_quantum_brute_force( # loop below therefore writes the weights into this uint64 buffer. out = np.empty_like(arrayf) - # Min weight of the block containing logical ops, with weights cast to uint64 as above. - weights = np.asarray(weight_func(arrayf[2**num_stabilizers :]), dtype=np.uint64) - min_weight = inplace_rowsum(weights).min(initial=num_bits) - if min_weight <= cutoff: - return int(min_weight) + # Until an operator is weighed, the lightest weight seen is one more than any operator can have, + # and an early exit must not return it. + min_weight = (num_bits if homogeneous else num_bits // 2) + 1 + + # Min weight of the block containing logical ops, if any, with weights cast to uint64 as above. + logical_block = arrayf[2**num_stabilizers :] + if len(logical_block): + weights = np.asarray(weight_func(logical_block), dtype=np.uint64) + min_weight = inplace_rowsum(weights).min() + if min_weight <= cutoff: + return int(min_weight) # Sweep over every remaining logical-op combination and, nested, every remaining stabilizer # combination, on top of the vectorized block above. The sweep uses a reflected Gray code: diff --git a/src/qldpc/codes/distance/_search/brute_force_test.py b/src/qldpc/codes/distance/_search/brute_force_test.py index facc4345f..ab493bc80 100644 --- a/src/qldpc/codes/distance/_search/brute_force_test.py +++ b/src/qldpc/codes/distance/_search/brute_force_test.py @@ -17,6 +17,7 @@ from ...conftest import Tune from .. import _binary as binary from ..conftest import count_set_bits +from . import brute_force @pytest.mark.parametrize("batch_exponent", [1, 3, 13]) @@ -169,3 +170,34 @@ def test_cutoff_early_exit(tune: Tune) -> None: method="brute_force", ) assert distance == 1 + + +@pytest.mark.parametrize("homogeneous", [True, False], ids=["hamming", "symplectic"]) +def test_early_exit_returns_only_weighed_operators(homogeneous: bool) -> None: + """An early exit returns the weight of an operator that the search has weighed. + + With a batch of one word, the vectorized block holds no logical operator, so the search has + weighed none before its sweep, and a cutoff at or above the number of columns must not end it + there. + """ + num_qubits = 3 + logical_op = [1, 0, 0] if homogeneous else [1, 0, 0, 0, 0, 0] + stabilizer = [0, 1, 1] if homogeneous else [0, 1, 1, 0, 0, 0] + distance = brute_force.get_distance_quantum_brute_force( + [logical_op], + [stabilizer], + cutoff=len(logical_op), + batch_exponent=0, + homogeneous=homogeneous, + ) + assert distance == 1 + + # Without logical operators, the result exceeds every weight, whatever the cutoff. + distance = brute_force.get_distance_quantum_brute_force( + np.empty((0, len(logical_op)), dtype=np.uint8), + [stabilizer], + cutoff=100, + batch_exponent=0, + homogeneous=homogeneous, + ) + assert distance == num_qubits + 1 diff --git a/src/qldpc/codes/distance/_search/completion_bounds.py b/src/qldpc/codes/distance/_search/completion_bounds.py index 28dd9256a..ce4d26373 100644 --- a/src/qldpc/codes/distance/_search/completion_bounds.py +++ b/src/qldpc/codes/distance/_search/completion_bounds.py @@ -18,7 +18,7 @@ for a minimum weight uses that bound until it finds a word within the maximum weight, and then, when each scan starts, the Hamming weight of its lightest word, or one less where every word of the coset has its parity. Once a scan finds a lighter word, the tables skip at the Hamming weight of the -lightest word that it has found (see ``coset_leader_search._search_coset_for_min_weight_by_tails``). +lightest word that it has found (see ``coset_leader_search._search_coset_for_min_weight``). The bases of the coset searches need not be information sets, so in place of ``w`` they bound the ones of a completion on pivot columns (see ``affine_search.scan_coset_layer_by_tails_with_tables``). diff --git a/src/qldpc/codes/distance/_search/connected_support.py b/src/qldpc/codes/distance/_search/connected_support.py index 0fc085278..c383289b9 100644 --- a/src/qldpc/codes/distance/_search/connected_support.py +++ b/src/qldpc/codes/distance/_search/connected_support.py @@ -105,8 +105,9 @@ def get_distance_connected_support( if not problems and upper_bound is None: raise ValueError("Require at least one problem or an upper bound") - # A problem without witnesses only bounds the result, and a known bound within the cutoff ends - # the search before it starts, as for get_distance_brouwer_zimmermann_many. + # A problem without witnesses only bounds the result, unlike in + # get_distance_brouwer_zimmermann_many, which needs a witness in every problem. A known bound + # within the cutoff ends the search before it starts, as it does there. tuning = get_tuning() if tuning is None else tuning automorphism_minimum_states = tuning.automorphism_minimum_states best_weight = upper_bound diff --git a/src/qldpc/codes/distance/_search/coset_leader_search.py b/src/qldpc/codes/distance/_search/coset_leader_search.py index 292876345..d73a70e3c 100644 --- a/src/qldpc/codes/distance/_search/coset_leader_search.py +++ b/src/qldpc/codes/distance/_search/coset_leader_search.py @@ -20,7 +20,6 @@ import numpy as np import numpy.typing as npt -from qldpc import math from qldpc._binary import ( get_kernel, pack_bits, @@ -54,6 +53,16 @@ # row that caches a stopping step (see ``_get_planning``). _PlanningBuffer = npt.NDArray[np.int64] +# The fields of a ``_CosetTails`` record, which the compiled searches take as the record itself or +# as a plain tuple. +_TailFields = tuple[ + npt.NDArray[np.uint64], + npt.NDArray[np.uint64], + npt.NDArray[np.uint8], + npt.NDArray[np.int64], + int, +] + #################################################################################################### # coset-leader minimum-weight bounds @@ -192,7 +201,7 @@ def __init__( self._support_searches: tuple[_ConnectedSupportCosetSearch, ...] = () self.css_weight = css_weight assert_binary(stabilizers, "stabilizers") - matrix = as_binary_matrix(stabilizers) + matrix = as_binary_matrix(stabilizers, "stabilizers") self.num_qubits = len(matrix) n = self.num_qubits if not n or matrix.shape[1] != 2 * n: @@ -200,7 +209,7 @@ def __init__( # The syndrome of a packed Pauli ``e`` is ``checks @ e``, where ``checks`` swaps the X and Z # halves of the stabilizers, and the stabilizers commute exactly when ``checks @ matrix.T`` # is zero. - checks = (math.symplectic_conjugate(matrix.astype(int)) % 2).astype(np.uint8) + checks = np.roll(matrix, n, axis=1) if len(row_reduce_binary(matrix)[2]) != n or np.any((checks @ matrix.T) % 2): raise ValueError("Stabilizers must be independent and commute") @@ -520,7 +529,7 @@ class _CosetSearchTables(NamedTuple): class _CosetTails(NamedTuple): - """The tails with which ``_search_coset_by_tails`` weighs the words of a family, and its tables. + """The tails with which ``_search_coset`` weighs the words of a family by tails, and its tables. The tail of an embedded word on a set of the family holds every column of the word but the pivot columns of the set, including columns on which no row of the set has a one. The tails are @@ -647,15 +656,10 @@ def get_min_weight_bound(self, syndrome: int, max_weight: int) -> MinWeightBound # Search the coset with the tables of the route of the maximum weight, by tails on a route # by tails, and keep the lightest word that the search weighs in ``best``. - tables, tails = self._routes.get_route(max_weight) + tables, tails = self._get_route(max_weight) best = np.empty(tables.inverse.shape[1], dtype=np.uint64) packed_syndrome = _int_to_words(syndrome, max(1, (self._num_syndrome_bits + 63) // 64)) - if tails is None: - lower_bound = _search_coset(packed_syndrome, max_weight, best, *tables, self.scale) - else: - lower_bound = _search_coset_by_tails( - packed_syndrome, max_weight, best, *tables, *tails, self.scale - ) + lower_bound = _search_coset(packed_syndrome, max_weight, best, *tables, tails, self.scale) # A search that finds a light enough word returns zero, but every word of the coset weighs # at least two, since the coset has no word of weight at most one. return MinWeightBound(lower_bound <= max_weight, self.unembed(best), max(2, lower_bound)) @@ -673,16 +677,11 @@ def find_first_not_within_max_weight( ``MinWeightBound`` for its coset, with the lightest word that the search weighed, or None if there is none. """ - tables, tails = self._routes.get_route(max_weight) + tables, tails = self._get_route(max_weight) best = np.empty(tables.inverse.shape[1], dtype=np.uint64) - if tails is None: - index, lower_bound = _find_first_not_within_max_weight( - syndromes, max_weight, best, *tables, self.scale - ) - else: - index, lower_bound = _find_first_not_within_max_weight_by_tails( - syndromes, max_weight, best, *tables, *tails, self.scale - ) + index, lower_bound = _find_first_not_within_max_weight( + syndromes, max_weight, best, *tables, tails, self.scale + ) if index == len(syndromes): return index, None return index, MinWeightBound(False, self.unembed(best), int(lower_bound)) @@ -705,12 +704,8 @@ def get_min_weights( that is larger (see ``_search_coset_for_min_weight``). """ max_weight = min(max_weight, self._num_qubits) - tables, tails = self._routes.get_route(max_weight) - if tails is None: - return _find_min_weights(syndromes, max_weight, floors, *tables, self.scale) - return _find_min_weights_by_tails( - syndromes, max_weight, floors, *tables, *tails, self.scale - ) + tables, tails = self._get_route(max_weight) + return _find_min_weights(syndromes, max_weight, floors, *tables, tails, self.scale) def unembed(self, embedded_word: npt.NDArray[np.uint64]) -> int: """Convert a packed embedded word into a word, without the columns that the embedding adds. @@ -720,6 +715,15 @@ def unembed(self, embedded_word: npt.NDArray[np.uint64]) -> int: """ return _packed_row_int(embedded_word) & ((1 << self._num_word_bits) - 1) + def _get_route(self, max_weight: int) -> tuple[_CosetSearchTables, _TailFields | None]: + """The tables of the route of a maximum weight, and the fields of its tails, if any. + + The fields are a plain tuple, which numba types at each call of a compiled search several + times faster than the named tuple ``_CosetTails``. + """ + tables, tails = self._routes.get_route(max_weight) + return tables, None if tails is None else tails[:] + @property def _num_word_bits(self) -> int: """The number of bits of a word: ``2 n`` for Paulis, and ``n`` for X or Z parts.""" @@ -1464,18 +1468,29 @@ def _find_first_not_within_max_weight( pivots: npt.NDArray[np.int64], planning: _PlanningBuffer, x_mask: npt.NDArray[np.uint64], + tails: _TailFields | None, scale: int, ) -> tuple[int, int]: """Search the coset of each syndrome in turn with ``_search_coset``. - The tables, from ``inverse`` to ``x_mask``, are those of ``_search_coset``. Return the index of - the first syndrome whose coset has no word of weight at most ``max_weight``, and the lower bound - that ``_search_coset`` returns for it, with its lightest word left in ``best``. If there is no - such syndrome, return the number of syndromes and zero. + The tables, from ``inverse`` to ``x_mask``, and ``tails`` are those of ``_search_coset``. + Return the index of the first syndrome whose coset has no word of weight at most + ``max_weight``, and the lower bound that ``_search_coset`` returns for it, with the lightest + word that it weighed left in ``best``. If there is no such syndrome, return the number of + syndromes and zero. """ for index in range(len(syndromes)): lower_bound = _search_coset( - syndromes[index], max_weight, best, inverse, generators, pivots, planning, x_mask, scale + syndromes[index], + max_weight, + best, + inverse, + generators, + pivots, + planning, + x_mask, + tails, + scale, ) if lower_bound > max_weight: return index, lower_bound @@ -1492,6 +1507,7 @@ def _search_coset( pivots: npt.NDArray[np.int64], planning: _PlanningBuffer, x_mask: npt.NDArray[np.uint64], + tails: _TailFields | None, scale: int, ) -> int: """Search the coset of a syndrome for a word of weight at most ``max_weight``. @@ -1509,7 +1525,11 @@ def _search_coset( the embedded word of a word whose syndrome has only bit ``j`` set. ``generators`` and ``pivots`` are the information sets of a family of the embedded generators of the code, with pivots padded with -1. ``planning`` is the planning buffer of the search with this family, - whose cached stopping steps a query may update. + whose cached stopping steps a query may update. ``tails`` is None, or, on a route by tails + (see ``coset_routes``), the fields of the ``_CosetTails`` of the family, and then ``scale`` is + one or two, since the search weighs words by their Hamming weight. Numba compiles the search + separately for None and for tails, and leaves out the branches that the type of ``tails`` rules + out, so choosing between them costs nothing for each word. The embedded words of the coset form a coset of the embedded code. The search starts from the sum of the rows of ``inverse`` for the set bits of the syndrome, a word with that syndrome. The @@ -1528,118 +1548,23 @@ def _search_coset( weighs at least the bound of step ``(w, j)`` in its table. So the search also stops at its stopping step, the first step whose bound exceeds ``max_weight``, which ``_get_stopping_step`` finds. It skips each set that adds nothing to the bound of that step (see - ``_search_from_syndrome_word``). + ``_search_from_syndrome_word``). A search by tails visits the same words in the same order, + except those that the tables of ``tails`` prove to lie beyond ``max_weight``, which it skips. + So it finds the same first word of weight at most ``max_weight``, if any. Returns: A lower bound on the weight of every word in the coset: zero, which is at most ``max_weight``, if ``best`` has weight at most ``max_weight``, and otherwise a bound that - exceeds ``max_weight``. + exceeds ``max_weight`` (see ``_search_from_syndrome_word``). """ best_weight, parity = _get_syndrome_word(packed_syndrome, best, inverse, x_mask, scale) if best_weight <= max_weight: return 0 table = _get_planning_table(planning, parity) stop_layer, stop_num_completed, stop_bound = _get_stopping_step(table, max_weight) - lower_bound, _ = _search_from_syndrome_word( - max_weight, - stop_layer, - stop_num_completed, - stop_bound, - best, - best_weight, - generators, - pivots, - x_mask, - scale, - ) - return lower_bound - - -@always_compiled -def _find_first_not_within_max_weight_by_tails( - syndromes: npt.NDArray[np.uint64], - max_weight: int, - best: npt.NDArray[np.uint64], - inverse: npt.NDArray[np.uint64], - generators: npt.NDArray[np.uint64], - pivots: npt.NDArray[np.int64], - planning: _PlanningBuffer, - x_mask: npt.NDArray[np.uint64], - tails: npt.NDArray[np.uint64], - inverse_tails: npt.NDArray[np.uint64], - completion_tables: npt.NDArray[np.uint8], - layout: npt.NDArray[np.int64], - block_bits: int, - scale: int, -) -> tuple[int, int]: - """Search the coset of each syndrome in turn with ``_search_coset_by_tails``. - - The arguments are those of ``_find_first_not_within_max_weight``, with the tails and tables of - ``_CosetTails``, and so is the result. Its lower bound and the word in ``best``, however, are - those that ``_search_coset_by_tails`` gives for the first syndrome without a word within - ``max_weight``. - """ - for index in range(len(syndromes)): - lower_bound = _search_coset_by_tails( - syndromes[index], - max_weight, - best, - inverse, - generators, - pivots, - planning, - x_mask, - tails, - inverse_tails, - completion_tables, - layout, - block_bits, - scale, - ) - if lower_bound > max_weight: - return index, lower_bound - return len(syndromes), 0 - - -@always_compiled(inline=True) -def _search_coset_by_tails( - packed_syndrome: npt.NDArray[np.uint64], - max_weight: int, - best: npt.NDArray[np.uint64], - inverse: npt.NDArray[np.uint64], - generators: npt.NDArray[np.uint64], - pivots: npt.NDArray[np.int64], - planning: _PlanningBuffer, - x_mask: npt.NDArray[np.uint64], - tails: npt.NDArray[np.uint64], - inverse_tails: npt.NDArray[np.uint64], - completion_tables: npt.NDArray[np.uint8], - layout: npt.NDArray[np.int64], - block_bits: int, - scale: int, -) -> int: - """Search the coset of a syndrome as ``_search_coset`` does, but by tails and with tables. - - The arguments are those of ``_search_coset``, with the tails and tables of ``_CosetTails`` of - the family of ``generators``, and ``scale`` is one or two, since the search weighs words by - their Hamming weight. The search visits the words of ``_search_coset`` in the same order, - except those that the tables prove to lie beyond ``max_weight``, which it skips (see - ``_search_by_tails_from_syndrome_word``). So it finds the same first word of weight at most - ``max_weight``, if any. Otherwise, it keeps the lightest word that it weighs in ``best``. - - Returns: - Zero if the search finds a word of weight at most ``max_weight``, and otherwise the lower - bound of ``_search_by_tails_from_syndrome_word``, which is the least weight beyond - ``max_weight`` that a word of the coset can have if the search skipped a word. - """ - best_weight, parity = _get_syndrome_word(packed_syndrome, best, inverse, x_mask, scale) - if best_weight <= max_weight: - return 0 - table = _get_planning_table(planning, parity) - stop_layer, stop_num_completed, stop_bound = _get_stopping_step(table, max_weight) - # The search rounds to the parity of the coset if its planning buffer has a table for each - # parity (see ``_get_planning_table``). - lower_bound, _, _ = _search_by_tails_from_syndrome_word( + # The search rounds to the parity of the coset if it weighs Hamming weights and its planning + # buffer has a table for each parity (see ``_get_planning_table``). + lower_bound, _, _ = _search_from_syndrome_word( packed_syndrome, max_weight, stop_layer, @@ -1649,13 +1574,10 @@ def _search_coset_by_tails( best_weight, generators, pivots, + x_mask, tails, - inverse_tails, - completion_tables, - layout, - block_bits, scale, - parity if len(planning) == 2 else -1, + parity if scale != 0 and len(planning) == 2 else -1, ) return lower_bound @@ -1684,10 +1606,10 @@ def _get_stopping_step(table: npt.NDArray[np.int64], max_weight: int) -> tuple[i ``[w, j - 1]``, and for a nonnegative maximum weight, step ``(w, j)`` is the stopping step exactly if its bound exceeds the maximum weight and entry ``[w, j - 1]`` does not. The step cached in the table (see ``_get_planning``) is reused if it passes this check. Otherwise, - ``information_sets.plan_steps`` finds the stopping step by binary search, and it is cached in - its place. A negative maximum weight fails the check, and the binary search returns the first - step. If no bound exceeds the maximum weight, which happens only when the maximum weight - reaches the largest native integer, the binary search returns the last step. + ``information_sets.find_stopping_step`` finds the stopping step by binary search, and it is + cached in its place. A negative maximum weight fails the check, and the binary search returns + the first step. If no bound exceeds the maximum weight, which happens only when the maximum + weight reaches the largest native integer, the binary search returns the last step. Oracles that share a cached search share its planning buffer. Compiled calls hold the global interpreter lock, since numba compiles them without ``nogil``, so calls interleave only when @@ -1705,11 +1627,8 @@ def _get_stopping_step(table: npt.NDArray[np.int64], max_weight: int) -> tuple[i max_weight < bounds[layer, num_completed] and bounds[layer, num_completed - 1] <= max_weight ) if not is_stopping_step: - # The search needs the step alone. Empty binomials ask ``plan_steps`` for the step only, - # and a plan that leaves out no sets reads no ranks, so both are empty placeholders. - no_values = np.empty(0, dtype=np.int64) - layer, num_completed, _ = information_sets.plan_steps( - bounds, max_weight, 0, len(bounds) - 1, no_values, no_values, False + layer, num_completed = information_sets.find_stopping_step( + bounds, max_weight, 0, len(bounds) - 1 ) table[-1, 0], table[-1, 1] = layer, num_completed return layer, num_completed, int(bounds[layer, num_completed]) @@ -1717,6 +1636,7 @@ def _get_stopping_step(table: npt.NDArray[np.int64], max_weight: int) -> tuple[i @always_compiled def _search_from_syndrome_word( + packed_syndrome: npt.NDArray[np.uint64], max_weight: int, stop_layer: int, stop_num_completed: int, @@ -1726,197 +1646,160 @@ def _search_from_syndrome_word( generators: npt.NDArray[np.uint64], pivots: npt.NDArray[np.int64], x_mask: npt.NDArray[np.uint64], + tails: _TailFields | None, scale: int, -) -> tuple[int, int]: + parity: int, +) -> tuple[int, int, int]: """Search a coset, as ``_search_coset`` does, from the word of its syndrome. - ``best`` holds that word, of weight ``best_weight`` above ``max_weight``. The stopping step and - its bound are those of ``_get_stopping_step``. A set of rank deficiency ``s`` adds to the bound - of a step only once it has completed layer ``s``, so a set that would not have completed it by - the stopping step adds nothing, and the search skips it. These are the unused sets of + ``best`` holds that word, of weight ``best_weight`` above ``max_weight``, and + ``packed_syndrome`` is its syndrome. The stopping step and its bound are those of + ``_get_stopping_step``. ``parity`` is the parity of the Hamming weights of the words of the + coset if the search weighs Hamming weights and rounds its bounds to that parity (see + ``_get_search_tables``), and otherwise -1. A set of rank deficiency ``s`` adds to the bound of + a step only once it has completed layer ``s``, so a set that would not have completed it by the + stopping step adds nothing, and the search skips it. These are the unused sets of ``information_sets.plan_steps``, whose work for this step, from layer zero, counts the words that the search weighs after the word of the syndrome if it finds no light enough word. - Returns: - The lower bound of ``_search_coset``, and the number of words weighed after the word of the - syndrome. Only tests read the number, to check which sets a search skips and that a - search by tails weighs as many words. - """ - num_sets, num_rows, width = generators.shape - offsets = _get_offsets(best, generators, pivots) - - # Visit each layer one information set at a time, up to the stopping step. The scans of layers - # take work space for one set, and weigh Hamming weights, which are ``scale`` times the weights, - # or CSS weights if ``scale`` is zero. - weighed = 0 - widths = np.full(1, width, dtype=np.int64) - support = np.empty(num_rows, dtype=np.int64) - prefixes = np.empty((num_rows, 1, width), dtype=np.uint64) - for layer in range(stop_layer + 1): - for set_index in range(num_sets if layer < stop_layer else stop_num_completed): - if information_sets.is_unused_set(pivots, set_index, stop_layer, stop_num_completed): - continue - if layer == 0: - best_weight = _keep_lighter(offsets[set_index], best, best_weight, x_mask, scale) - weighed += 1 - elif scale == 0: - best_weight, count = affine_search.scan_coset_layer_by_css_weight( - offsets[set_index], - generators[set_index : set_index + 1], - widths, - support[:layer], - prefixes, - best_weight, - max_weight, - x_mask, - best, - ) - weighed += count - else: - best_ones, count = affine_search.scan_coset_layer( - offsets[set_index], - generators[set_index : set_index + 1], - widths, - support[:layer], - prefixes, - scale * best_weight, - scale * max_weight, - best, - ) - best_weight = best_ones // scale - weighed += count - if best_weight <= max_weight: - return 0, weighed - return min(best_weight, stop_bound), weighed - - -@always_compiled -def _search_by_tails_from_syndrome_word( - packed_syndrome: npt.NDArray[np.uint64], - max_weight: int, - stop_layer: int, - stop_num_completed: int, - stop_bound: int, - best: npt.NDArray[np.uint64], - best_weight: int, - generators: npt.NDArray[np.uint64], - pivots: npt.NDArray[np.int64], - tails: npt.NDArray[np.uint64], - inverse_tails: npt.NDArray[np.uint64], - completion_tables: npt.NDArray[np.uint8], - layout: npt.NDArray[np.int64], - block_bits: int, - scale: int, - parity: int, -) -> tuple[int, int, int]: - """Search a coset as ``_search_from_syndrome_word`` does, but by tails and with tables. - - The arguments are those of ``_search_from_syndrome_word``, with ``packed_syndrome``, the - syndrome of the word in ``best``, and the tails and tables of ``_CosetTails`` rather than - ``x_mask``. ``scale`` is one or two. ``parity`` is the parity of the Hamming weights of the - words of the coset if the search rounds its bounds to it (see ``_get_search_tables``), and - otherwise -1. The search visits the same words in the same order, except those that the tables - skip, and keeps the first of the lightest that it weighs, but it weighs a word by its rows below - the rank of its set and its tail. It keeps the set, the layer, and the support of the lightest - word, and writes that word into ``best`` before it returns. - - In each layer of two or more rows of a set with tables, the search skips the completions of a - prefix whose bound reaches the skip bound: the least bound that proves a word of the coset - heavier than ``max_weight``, one more than the greatest Hamming weight at most ``scale * - max_weight`` that a word of the coset can have (see + A search by tails (see ``_search_coset``) weighs a word by its rows below the rank of its set + and its tail. It keeps the set, the layer, and the support of the lightest word, and writes + that word into ``best`` before it returns. In each layer of two or more rows of a set with + tables, it skips the completions of a prefix whose bound reaches the skip bound: the least bound + that proves a word of the coset heavier than ``max_weight``, one more than the greatest Hamming + weight at most ``scale * max_weight`` that a word of the coset can have (see ``affine_search.scan_coset_layer_by_tails_with_tables``). So it finds the same first word of weight at most ``max_weight``, if any. If it finds none and skipped any word, it bounds the coset by the least weight beyond ``max_weight`` that a word of the coset can have, which bounds the skipped words, the words that it weighed, and the words beyond the stopping step. Returns: - The lower bound of ``_search_coset`` if the search finds a word of weight at most - ``max_weight`` or skips no word, and otherwise the least weight beyond ``max_weight`` that a - word of the coset can have; and the numbers of words weighed after the word of the syndrome - and of supports skipped. Only tests read the numbers. + A lower bound on the weight of every word in the coset: zero if the search finds a word of + weight at most ``max_weight``, the least weight beyond ``max_weight`` that a word of the + coset can have if it skipped any word, and otherwise the lesser of the weight of the + lightest word that it weighed and the bound of the stopping step. Also the numbers of words + weighed after the word of the syndrome and of supports skipped. Only tests read the + numbers, to check which sets a search skips and that a search by tails weighs as many words. """ - num_sets, num_rows, _ = generators.shape - tail_width = tails.shape[2] - ranks, offset_tails = _get_offset_tails(packed_syndrome, best, pivots, tails, inverse_tails) + # Reduce the word of the syndrome to the offset of each information set, or, by tails, find the + # rank of each set and the tail of its offset. + num_sets, num_rows, width = generators.shape + if tails is None: + offsets = _get_offsets(best, generators, pivots) + else: + tail_rows, inverse_tails, completion_tables, layout, block_bits = tails + ranks, offsets = _get_offset_tails(packed_syndrome, best, pivots, tail_rows, inverse_tails) + width = tail_rows.shape[2] - # The scans of layers weigh Hamming weights, which are ``scale`` times the weights. - best_ones, cutoff = scale * best_weight, scale * max_weight + # The scans weigh ``factor`` units for each unit of weight: Hamming weights, which are ``scale`` + # times the weights, or CSS weights if ``scale`` is zero. + factor = max(1, scale) + best_units, cutoff = factor * best_weight, factor * max_weight skip_bound = _get_skip_bound(cutoff, parity) - # Visit each layer one information set at a time, up to the stopping step, and skip the sets - # that ``_search_from_syndrome_word`` skips. + # Visit each layer one information set at a time, up to the stopping step. The scans take work + # space for one set. Those by tails keep the support of each lighter word in ``best_support``, + # and the others keep the word itself in ``best``. best_set, best_layer = -1, 0 weighed, skipped = 0, 0 - widths = np.full(1, tail_width, dtype=np.int64) + widths = np.full(1, width, dtype=np.int64) support = np.empty(num_rows, dtype=np.int64) best_support = np.empty(num_rows, dtype=np.int64) - prefixes = np.empty((num_rows, 1, tail_width), dtype=np.uint64) + prefixes = np.empty((num_rows, 1, width), dtype=np.uint64) for layer in range(stop_layer + 1): for set_index in range(num_sets if layer < stop_layer else stop_num_completed): if information_sets.is_unused_set(pivots, set_index, stop_layer, stop_num_completed): continue + set_skipped = 0 if layer == 0: - ones = 0 - for index in range(tail_width): - ones += int(popcount_word(offset_tails[set_index, index])) - count = 1 + # The offset is zero on the pivot columns, so its tail has its Hamming weight. + weight, ones = _get_weight(offsets[set_index], x_mask, scale) + found, count = weight if scale == 0 else ones, 1 + if tails is None and found < best_units: + best[:] = offsets[set_index] + elif tails is None: + if scale == 0: + found, count = affine_search.scan_coset_layer_by_css_weight( + offsets[set_index], + generators[set_index : set_index + 1], + widths, + support[:layer], + prefixes, + best_units, + cutoff, + x_mask, + best, + ) + else: + found, count = affine_search.scan_coset_layer( + offsets[set_index], + generators[set_index : set_index + 1], + widths, + support[:layer], + prefixes, + best_units, + cutoff, + best, + ) elif layer >= 2 and layout[set_index, 1]: # The tables of the set start at ``layout[set_index, 0]``, with a row of tables for # each row of the set but the last. num_blocks = layout[set_index, 1] block_offsets = layout[set_index, 2 : num_blocks + 3] start = layout[set_index, 0] - ones, count, set_skipped = affine_search.scan_coset_layer_by_tails_with_tables( - offset_tails[set_index], - tails[set_index : set_index + 1], + end = start + (num_rows - 1) * block_offsets[-1] + found, count, set_skipped = affine_search.scan_coset_layer_by_tails_with_tables( + offsets[set_index], + tail_rows[set_index : set_index + 1], ranks[set_index : set_index + 1], widths, support[:layer], prefixes, - best_ones, + best_units, cutoff, skip_bound, - completion_tables[start : start + (num_rows - 1) * block_offsets[-1]], + completion_tables[start:end], block_offsets, block_bits, best_support, ) - skipped += set_skipped else: - ones, count = affine_search.scan_coset_layer_by_tails( - offset_tails[set_index], - tails[set_index : set_index + 1], + found, count = affine_search.scan_coset_layer_by_tails( + offsets[set_index], + tail_rows[set_index : set_index + 1], ranks[set_index : set_index + 1], widths, support[:layer], prefixes, - best_ones, + best_units, cutoff, best_support, ) weighed += count - if ones < best_ones: - best_ones, best_set, best_layer = ones, set_index, layer - if best_ones <= cutoff: + skipped += set_skipped + if found < best_units: + best_units, best_set, best_layer = found, set_index, layer + if best_units <= cutoff: break - if best_ones <= cutoff: + if best_units <= cutoff: break - # Write the lightest word into ``best``, which holds the word of the syndrome until then. - if best_set >= 0: + # By tails, write the lightest word into ``best``, which holds the word of the syndrome until + # then. + if tails is not None and best_set >= 0: _rebuild_word(best, generators[best_set], pivots[best_set], best_support[:best_layer]) # A search that finds no light enough word but skipped words bounds the coset by the least # weight beyond the maximum weight that a word of the coset can have: the skip bound rounded up # to the parity of the coset, then divided by the scale and rounded up, as the bounds of the # steps are. - if best_ones <= cutoff: + if best_units <= cutoff: lower_bound = 0 elif skipped: rounded = skip_bound + (skip_bound - parity) % 2 if parity >= 0 else skip_bound - lower_bound = -(-rounded // scale) + lower_bound = -(-rounded // factor) else: - lower_bound = min(best_ones // scale, stop_bound) + lower_bound = min(best_units // factor, stop_bound) return lower_bound, weighed, skipped @@ -1930,6 +1813,7 @@ def _find_min_weights( pivots: npt.NDArray[np.int64], planning: _PlanningBuffer, x_mask: npt.NDArray[np.uint64], + tails: _TailFields | None, scale: int, ) -> npt.NDArray[np.int64]: """Search the coset of each syndrome in turn with ``_search_coset_for_min_weight``. @@ -1941,7 +1825,7 @@ def _find_min_weights( weights = np.empty(len(syndromes), dtype=np.int64) best = np.empty(inverse.shape[1], dtype=np.uint64) for index in range(len(syndromes)): - weight, _ = _search_coset_for_min_weight( + weight, _, _ = _search_coset_for_min_weight( syndromes[index], max_weight, floors[index], @@ -1951,6 +1835,7 @@ def _find_min_weights( pivots, planning, x_mask, + tails, scale, ) weights[index] = weight @@ -1968,8 +1853,9 @@ def _search_coset_for_min_weight( pivots: npt.NDArray[np.int64], planning: _PlanningBuffer, x_mask: npt.NDArray[np.uint64], + tails: _TailFields | None, scale: int, -) -> tuple[int, int]: +) -> tuple[int, int, int]: """Search the coset of a syndrome for its minimum weight, if that is at most ``max_weight``. The arguments are those of ``_search_coset``, with a floor of at most ``max_weight``, and @@ -1994,22 +1880,40 @@ def _search_coset_for_min_weight( best word, no word is one lighter, so the incumbent is one less and the scans keep the same words. This holds in a search of Pauli errors by physical weight, whose embedded words have even Hamming weights, and in a search of X or Z parts that rounds to the parity of each coset. - The lower incumbent lets the tables of ``_search_coset_for_min_weight_by_tails`` skip more. + The lower incumbent lets the tables of a search by tails skip more. + + A search by tails visits the same words in the same order, but weighs a word by its rows below + the rank of its set and its tail, as ``_search_from_syndrome_word`` does. In each layer of two + or more rows of a set with tables, it skips the completions of a prefix whose bound reaches the + best weight of the scan, which starts at its incumbent: the scan takes the skip bound of + ``max_weight``, which is at least the incumbent, as its skip bound (see + ``affine_search.scan_coset_layer_by_tails_with_tables``). Until the search finds a word within + ``max_weight``, the incumbent is that skip bound, at which a query by tails skips too. Then it + is the Hamming weight of the best word, or one less where every word of the coset has its + parity, so the search skips no word lighter than the best word. So the search gives the weight + of a search of whole words. For a coset without a word within ``max_weight``, it weighs and + skips the words that ``_search_from_syndrome_word`` weighs and skips at ``max_weight``. Returns: The best weight when the search stops, or ``max_weight + 1`` if it found no word within - ``max_weight``, and the number of words weighed after the word of the syndrome. The best - weight is the minimum weight of the coset, unless it is at most ``floor``. Only tests read - the number. + ``max_weight``, and the numbers of words weighed after the word of the syndrome and of + supports skipped. The best weight is the minimum weight of the coset, unless it is at most + ``floor``. Only tests read the numbers. """ - # Start from the word of the syndrome, which ends the search if it is within the floor, and - # reduce it to the offset of each information set. + # Start from the word of the syndrome, which ends the search if it is within the floor. Reduce + # it to the offset of each information set, or, by tails, find the rank of each set and the + # tail of its offset. best_weight, parity = _get_syndrome_word(packed_syndrome, best, inverse, x_mask, scale) if best_weight <= floor: - return best_weight, 0 + return best_weight, 0, 0 bounds = _get_planning_table(planning, parity)[:-1] num_sets, num_rows, width = generators.shape - offsets = _get_offsets(best, generators, pivots) + if tails is None: + offsets = _get_offsets(best, generators, pivots) + else: + tail_rows, inverse_tails, completion_tables, layout, block_bits = tails + ranks, offsets = _get_offset_tails(packed_syndrome, best, pivots, tail_rows, inverse_tails) + width = tail_rows.shape[2] # Find the skip bound and the parity slack of the incumbent, in the units of the scans, which # weigh ``factor`` units for each unit of weight. Only Hamming weights have a parity. Find @@ -2026,10 +1930,13 @@ def _search_coset_for_min_weight( # Visit each layer one information set at a time, up to the stopping step. The scans take work # space for one set. The search reaches a stopping step by the last layer, whose bounds are # ``NO_WORD_BOUND``, since its upper bound is at most ``max_weight + 1``, which - # ``_CosetSearch.get_min_weights`` keeps at most ``n + 1``. - weighed = 0 + # ``_CosetSearch.get_min_weights`` keeps at most ``n + 1``. The scans keep each lighter word in + # ``best``, or its support in ``best_support`` by tails, which are only work space, since the + # search returns no word. + weighed, skipped = 0, 0 widths = np.full(1, width, dtype=np.int64) support = np.empty(num_rows, dtype=np.int64) + best_support = np.empty(num_rows, dtype=np.int64) prefixes = np.empty((num_rows, 1, width), dtype=np.uint64) layer = 0 while True: @@ -2045,206 +1952,61 @@ def _search_coset_for_min_weight( incumbent = min(best_units - slack, skip_bound) step_floor = max(floor, bounds[layer, set_index]) cutoff = min(factor * step_floor, incumbent - 1) + set_skipped = 0 if layer == 0: + # The offset is zero on the pivot columns, so its tail has its Hamming weight. weight, ones = _get_weight(offsets[set_index], x_mask, scale) found, count = weight if scale == 0 else ones, 1 - elif scale == 0: - found, count = affine_search.scan_coset_layer_by_css_weight( - offsets[set_index], - generators[set_index : set_index + 1], - widths, - support[:layer], - prefixes, - incumbent, - cutoff, - x_mask, - best, - ) - else: - found, count = affine_search.scan_coset_layer( - offsets[set_index], - generators[set_index : set_index + 1], - widths, - support[:layer], - prefixes, - incumbent, - cutoff, - best, - ) - weighed += count - if found < incumbent: - best_units, best_weight = found, found // factor - if best_weight <= step_floor: - return best_weight, weighed - - # Stop at the stopping step, and move that step earlier if the best weight fell. - new_upper_bound = min(best_weight, max_weight + 1) - if bounds[layer, set_index + 1] >= new_upper_bound: - return new_upper_bound, weighed - if new_upper_bound < upper_bound: - upper_bound = new_upper_bound - stop_layer, stop_num_completed = _get_first_step_reaching(bounds, upper_bound) - layer += 1 - - -@always_compiled -def _find_min_weights_by_tails( - syndromes: npt.NDArray[np.uint64], - max_weight: int, - floors: npt.NDArray[np.int64], - inverse: npt.NDArray[np.uint64], - generators: npt.NDArray[np.uint64], - pivots: npt.NDArray[np.int64], - planning: _PlanningBuffer, - x_mask: npt.NDArray[np.uint64], - tails: npt.NDArray[np.uint64], - inverse_tails: npt.NDArray[np.uint64], - completion_tables: npt.NDArray[np.uint8], - layout: npt.NDArray[np.int64], - block_bits: int, - scale: int, -) -> npt.NDArray[np.int64]: - """Search the coset of each syndrome in turn with ``_search_coset_for_min_weight_by_tails``. - - The arguments are those of ``_find_min_weights``, with the tails and tables of ``_CosetTails``, - and so is the result. - """ - weights = np.empty(len(syndromes), dtype=np.int64) - best = np.empty(inverse.shape[1], dtype=np.uint64) - for index in range(len(syndromes)): - weight, _, _ = _search_coset_for_min_weight_by_tails( - syndromes[index], - max_weight, - floors[index], - best, - inverse, - generators, - pivots, - planning, - x_mask, - tails, - inverse_tails, - completion_tables, - layout, - block_bits, - scale, - ) - weights[index] = weight - return weights - - -@always_compiled(inline=True) -def _search_coset_for_min_weight_by_tails( - packed_syndrome: npt.NDArray[np.uint64], - max_weight: int, - floor: int, - best: npt.NDArray[np.uint64], - inverse: npt.NDArray[np.uint64], - generators: npt.NDArray[np.uint64], - pivots: npt.NDArray[np.int64], - planning: _PlanningBuffer, - x_mask: npt.NDArray[np.uint64], - tails: npt.NDArray[np.uint64], - inverse_tails: npt.NDArray[np.uint64], - completion_tables: npt.NDArray[np.uint8], - layout: npt.NDArray[np.int64], - block_bits: int, - scale: int, -) -> tuple[int, int, int]: - """Search the coset of a syndrome as ``_search_coset_for_min_weight`` does, but by tails. - - The arguments are those of ``_search_coset_for_min_weight``, with the tails and tables of - ``_CosetTails``, and ``scale`` is one or two. The search visits the words of - ``_search_coset_for_min_weight`` in the same order, but weighs a word by its rows below the rank - of its set and its tail, as ``_search_by_tails_from_syndrome_word`` does. - - In each layer of two or more rows of a set with tables, it skips the completions of a prefix - whose bound reaches the best weight of the scan, which starts at its incumbent (see - ``affine_search.scan_coset_layer_by_tails_with_tables_below_best``). Until it finds a word - within ``max_weight``, the incumbent is the skip bound of ``max_weight``, at which - ``_search_by_tails_from_syndrome_word`` skips too. Then it is the Hamming weight of the best - word, or one less where every word of the coset has its parity, so the search skips no word - lighter than the best word. So the search gives the weight of ``_search_coset_for_min_weight``. - For a coset without a word within ``max_weight``, it weighs and skips the words that - ``_search_by_tails_from_syndrome_word`` weighs and skips at ``max_weight``. - - Returns: - The weight of ``_search_coset_for_min_weight``, and the numbers of words weighed after the - word of the syndrome and of supports skipped. Only tests read the numbers. - """ - # Start from the word of the syndrome, which ends the search if it is within the floor, and - # find the tail of its offset on each information set. - best_weight, parity = _get_syndrome_word(packed_syndrome, best, inverse, x_mask, scale) - if best_weight <= floor: - return best_weight, 0, 0 - bounds = _get_planning_table(planning, parity)[:-1] - num_sets, num_rows, _ = generators.shape - tail_width = tails.shape[2] - ranks, offset_tails = _get_offset_tails(packed_syndrome, best, pivots, tails, inverse_tails) - - # Find the skip bound and the parity slack of the incumbent, in Hamming weights. Find also the - # upper bound ``U`` of the search, in units of weight, and its first stopping step (see - # ``_search_coset_for_min_weight``). - rounds = len(planning) == 2 - skip_bound = _get_skip_bound(scale * max_weight, parity if rounds else -1) - slack = 1 if scale == 2 or rounds else 0 - best_ones = scale * best_weight - upper_bound = min(best_weight, max_weight + 1) - stop_layer, stop_num_completed = _get_first_step_reaching(bounds, upper_bound) - - # Visit each layer one information set at a time, up to the stopping step, which the search - # reaches as ``_search_coset_for_min_weight`` does. The scans keep the support of each lighter - # word in ``best_support``, which is only work space, since the search returns no word. - weighed, skipped = 0, 0 - widths = np.full(1, tail_width, dtype=np.int64) - support = np.empty(num_rows, dtype=np.int64) - best_support = np.empty(num_rows, dtype=np.int64) - prefixes = np.empty((num_rows, 1, tail_width), dtype=np.uint64) - layer = 0 - while True: - for set_index in range(num_sets): - # Weigh the words of this step that are lighter than the incumbent, unless the set adds - # nothing to the bound of the stopping step, and return once the best weight is at most - # the floor or the bound of the step before this one. A scan returns once its best - # weight, which starts at the incumbent, is at most its cutoff, so the cutoff stays - # below the incumbent. - if not information_sets.is_unused_set( - pivots, set_index, stop_layer, stop_num_completed - ): - incumbent = min(best_ones - slack, skip_bound) - step_floor = max(floor, bounds[layer, set_index]) - cutoff = min(scale * step_floor, incumbent - 1) - if layer == 0: - found, count = 0, 1 - for index in range(tail_width): - found += int(popcount_word(offset_tails[set_index, index])) - elif layer >= 2 and layout[set_index, 1]: - # The tables of the set start at ``layout[set_index, 0]`` (see - # ``_search_by_tails_from_syndrome_word``). - num_blocks = layout[set_index, 1] - block_offsets = layout[set_index, 2 : num_blocks + 3] - start = layout[set_index, 0] - found, count, set_skipped = ( - affine_search.scan_coset_layer_by_tails_with_tables_below_best( - offset_tails[set_index], - tails[set_index : set_index + 1], - ranks[set_index : set_index + 1], + elif tails is None: + if scale == 0: + found, count = affine_search.scan_coset_layer_by_css_weight( + offsets[set_index], + generators[set_index : set_index + 1], widths, support[:layer], prefixes, incumbent, cutoff, - completion_tables[start : start + (num_rows - 1) * block_offsets[-1]], - block_offsets, - block_bits, - best_support, + x_mask, + best, ) + else: + found, count = affine_search.scan_coset_layer( + offsets[set_index], + generators[set_index : set_index + 1], + widths, + support[:layer], + prefixes, + incumbent, + cutoff, + best, + ) + elif layer >= 2 and layout[set_index, 1]: + # The tables of the set start at ``layout[set_index, 0]`` (see + # ``_search_from_syndrome_word``). + num_blocks = layout[set_index, 1] + block_offsets = layout[set_index, 2 : num_blocks + 3] + start = layout[set_index, 0] + end = start + (num_rows - 1) * block_offsets[-1] + found, count, set_skipped = affine_search.scan_coset_layer_by_tails_with_tables( + offsets[set_index], + tail_rows[set_index : set_index + 1], + ranks[set_index : set_index + 1], + widths, + support[:layer], + prefixes, + incumbent, + cutoff, + skip_bound, + completion_tables[start:end], + block_offsets, + block_bits, + best_support, ) - skipped += set_skipped else: found, count = affine_search.scan_coset_layer_by_tails( - offset_tails[set_index], - tails[set_index : set_index + 1], + offsets[set_index], + tail_rows[set_index : set_index + 1], ranks[set_index : set_index + 1], widths, support[:layer], @@ -2254,8 +2016,9 @@ def _search_coset_for_min_weight_by_tails( best_support, ) weighed += count + skipped += set_skipped if found < incumbent: - best_ones, best_weight = found, found // scale + best_units, best_weight = found, found // factor if best_weight <= step_floor: return best_weight, weighed, skipped @@ -2273,16 +2036,11 @@ def _search_coset_for_min_weight_by_tails( def _get_first_step_reaching(bounds: npt.NDArray[np.int64], target: int) -> tuple[int, int]: """The first step whose bound in a table of bounds of steps reaches ``target``. - It is the stopping step of ``information_sets.plan_steps`` for the limit ``target - 1``. Given - empty ranks and binomials, ``plan_steps`` returns only that step, with zero work. Unlike - ``_get_stopping_step``, this function caches no step, so the cached steps of the planning buffer - serve only the threshold queries of ``_search_coset`` and ``_search_coset_by_tails``. + It is the stopping step of ``information_sets.find_stopping_step`` for the limit + ``target - 1``. Unlike ``_get_stopping_step``, this function caches no step, so the cached + steps of the planning buffer serve only the queries of ``_search_coset``. """ - no_values = np.empty(0, dtype=np.int64) - layer, num_completed, _ = information_sets.plan_steps( - bounds, target - 1, 0, len(bounds) - 1, no_values, no_values, False - ) - return layer, num_completed + return information_sets.find_stopping_step(bounds, target - 1, 0, len(bounds) - 1) @always_compiled(inline=True) @@ -2311,11 +2069,11 @@ def _get_offset_tails( ) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.uint64]]: """Find the rank of each information set and the tail of its offset. - ``word`` is the word of the syndrome ``packed_syndrome``, and ``pivots``, ``tails``, and - ``inverse_tails`` are those of ``_search_coset_by_tails``. The rows of a set with pivots on - which the word has ones reduce the word to the offset, since every row is zero on the pivots of - the other rows. The tail of the word sums the tails of the rows of the inverse for the set bits - of the syndrome. + ``word`` is the word of the syndrome ``packed_syndrome``, ``pivots`` are those of + ``_search_coset``, and ``tails`` and ``inverse_tails`` are those of ``_CosetTails``. The rows + of a set with pivots on which the word has ones reduce the word to the offset, since every row + is zero on the pivots of the other rows. The tail of the word sums the tails of the rows of the + inverse for the set bits of the syndrome. """ num_sets, num_rows, tail_width = tails.shape ranks = np.zeros(num_sets, dtype=np.int64) @@ -2420,22 +2178,6 @@ def _xor_selected_rows( bits &= bits - np.uint64(1) -@always_compiled -def _keep_lighter( - word: npt.NDArray[np.uint64], - best: npt.NDArray[np.uint64], - best_weight: int, - x_mask: npt.NDArray[np.uint64], - scale: int, -) -> int: - """Keep ``word`` in ``best`` if it is lighter, and return the weight of ``best``.""" - weight, _ = _get_weight(word, x_mask, scale) - if weight >= best_weight: - return best_weight - best[:] = word - return weight - - @always_compiled def _get_weight( word: npt.NDArray[np.uint64], x_mask: npt.NDArray[np.uint64], scale: int diff --git a/src/qldpc/codes/distance/_search/coset_leader_search_test.py b/src/qldpc/codes/distance/_search/coset_leader_search_test.py index 23c4d524d..ae88bed08 100644 --- a/src/qldpc/codes/distance/_search/coset_leader_search_test.py +++ b/src/qldpc/codes/distance/_search/coset_leader_search_test.py @@ -786,8 +786,8 @@ def test_oracles_build_other_families_with_the_tuning_they_are_built_with( signature = inspect.signature(get_overlapping_family) def record(*args: Any, **kwargs: Any) -> information_sets.BasisFamily: - arguments = signature.bind(*args, **kwargs).arguments - calls.append((arguments["max_sets"], arguments["attempts_per_set"])) + policy = signature.bind(*args, **kwargs).arguments["policy"] + calls.append((policy.max_sets, policy.attempts_per_set)) return get_overlapping_family(*args, **kwargs) def record_tables( @@ -981,13 +981,13 @@ def test_cached_stopping_steps_are_reused_exactly_while_they_stay_first( zero, two, and four in a code of dimension four, for both weight factors. """ plans: list[int] = [] - plan_steps = information_sets.plan_steps + find_stopping_step = information_sets.find_stopping_step - def record(*args: Any) -> tuple[int, int, int]: + def record(*args: Any) -> tuple[int, int]: plans.append(args[1]) - return plan_steps(*args) + return find_stopping_step(*args) - monkeypatch.setattr(information_sets, "plan_steps", record) + monkeypatch.setattr(information_sets, "find_stopping_step", record) # Take the Python version even when the function is compiled, since only it calls the recording # function. function = coset_leader_search._get_stopping_step @@ -1287,8 +1287,9 @@ def test_searches_visit_only_the_sets_that_add_to_their_bounds(pytestconfig: pyt range(1 << len(checks)), range(num_columns + 1) ): best = np.empty(tables.inverse.shape[1], dtype=np.uint64) + packed_syndrome = _pack_syndromes([syndrome], len(checks))[0] best_weight, parity = coset_leader_search._get_syndrome_word( - _pack_syndromes([syndrome], len(checks))[0], best, tables.inverse, tables.x_mask, 1 + packed_syndrome, best, tables.inverse, tables.x_mask, 1 ) if best_weight <= max_weight: continue @@ -1296,7 +1297,8 @@ def test_searches_visit_only_the_sets_that_add_to_their_bounds(pytestconfig: pyt stop_layer, stop_num_completed, stop_bound = coset_leader_search._get_stopping_step( table, max_weight ) - lower_bound, weighed = coset_leader_search._search_from_syndrome_word( + lower_bound, weighed, _ = coset_leader_search._search_from_syndrome_word( + packed_syndrome, max_weight, stop_layer, stop_num_completed, @@ -1306,7 +1308,9 @@ def test_searches_visit_only_the_sets_that_add_to_their_bounds(pytestconfig: pyt tables.generators, tables.pivots, tables.x_mask, + None, 1, + parity if len(tables.planning) == 2 else -1, ) representative = unpack_bits(best[None], num_columns)[0].astype(int) assert _classical_syndrome(checks, representative) == syndrome @@ -1414,7 +1418,8 @@ def _assert_searches_by_tails_match_searches_of_words( stop_layer, stop_num_completed, stop_bound = coset_leader_search._get_stopping_step( coset_leader_search._get_planning_table(tables.planning, parity), max_weight ) - answer = coset_leader_search._search_from_syndrome_word( + *answer, _ = coset_leader_search._search_from_syndrome_word( + packed_syndrome, max_weight, stop_layer, stop_num_completed, @@ -1424,9 +1429,11 @@ def _assert_searches_by_tails_match_searches_of_words( tables.generators, tables.pivots, tables.x_mask, + None, scale, + parity if len(tables.planning) == 2 else -1, ) - *answer_by_tails, skipped = coset_leader_search._search_by_tails_from_syndrome_word( + *answer_by_tails, skipped = coset_leader_search._search_from_syndrome_word( packed_syndrome, max_weight, stop_layer, @@ -1436,7 +1443,8 @@ def _assert_searches_by_tails_match_searches_of_words( best_weight, tables.generators, tables.pivots, - *tails, + tables.x_mask, + tails, scale, parity if len(tables.planning) == 2 else -1, ) @@ -1475,10 +1483,10 @@ def _assert_searches_by_tails_match_searches_of_words( while True: words = np.empty((2, tables.inverse.shape[1]), dtype=np.uint64) found = coset_leader_search._find_first_not_within_max_weight( - packed_syndromes[start:], max_weight, words[0], *tables, scale + packed_syndromes[start:], max_weight, words[0], *tables, None, scale ) - found_by_tails = coset_leader_search._find_first_not_within_max_weight_by_tails( - packed_syndromes[start:], max_weight, words[1], *tables, *tails, scale + found_by_tails = coset_leader_search._find_first_not_within_max_weight( + packed_syndromes[start:], max_weight, words[1], *tables, tails, scale ) assert found_by_tails[0] == found[0] if start + found[0] == len(packed_syndromes): @@ -1975,41 +1983,25 @@ def test_min_weight_searches_above_max_weight_examine_as_queries( stop_layer, stop_num_completed, stop_bound = coset_leader_search._get_stopping_step( coset_leader_search._get_planning_table(tables.planning, parity), max_weight ) - if tails is None: - _, *expected = coset_leader_search._search_from_syndrome_word( - max_weight, - stop_layer, - stop_num_completed, - stop_bound, - best, - best_weight, - tables.generators, - tables.pivots, - tables.x_mask, - search.scale, - ) - _, *found = coset_leader_search._search_coset_for_min_weight( - packed_syndrome, max_weight, 0, best, *tables, search.scale - ) - else: - _, *expected = coset_leader_search._search_by_tails_from_syndrome_word( - packed_syndrome, - max_weight, - stop_layer, - stop_num_completed, - stop_bound, - best, - best_weight, - tables.generators, - tables.pivots, - *tails, - search.scale, - parity if len(tables.planning) == 2 else -1, - ) - _, *found = coset_leader_search._search_coset_for_min_weight_by_tails( - packed_syndrome, max_weight, 0, best, *tables, *tails, search.scale - ) - skipped += expected[1] + _, *expected = coset_leader_search._search_from_syndrome_word( + packed_syndrome, + max_weight, + stop_layer, + stop_num_completed, + stop_bound, + best, + best_weight, + tables.generators, + tables.pivots, + tables.x_mask, + tails, + search.scale, + parity if len(tables.planning) == 2 else -1, + ) + _, *found = coset_leader_search._search_coset_for_min_weight( + packed_syndrome, max_weight, 0, best, *tables, tails, search.scale + ) + skipped += expected[1] assert found == expected assert skipped or not (tuning and must_skip) diff --git a/src/qldpc/codes/distance/_search/coset_routes.py b/src/qldpc/codes/distance/_search/coset_routes.py index 23b83bd56..fe1c83b95 100644 --- a/src/qldpc/codes/distance/_search/coset_routes.py +++ b/src/qldpc/codes/distance/_search/coset_routes.py @@ -294,13 +294,10 @@ def _get_sets_needing_tables( tables if the predicate considers its last layer. """ num_rows = len(self._basis) - # The routes need the stopping step alone. Empty binomials ask ``plan_steps`` for the step - # only, and a plan that leaves out no sets reads no ranks, so both are empty placeholders. - no_values = np.empty(0, dtype=np.int64) sets = set() for table in tables.planning: - stop_layer, stop_num_completed, _ = information_sets.plan_steps( - table[:-1], max_weight, 0, len(table) - 2, no_values, no_values, False + stop_layer, stop_num_completed = information_sets.find_stopping_step( + table[:-1], max_weight, 0, len(table) - 2 ) for set_index in range(len(instance.ranks)): last_layer = information_sets.get_completed_layer( @@ -333,7 +330,7 @@ def _get_candidates( np.empty((len(self._basis), 0), dtype=np.uint8), maximum if len(maximum.ranks) else disjoint, weight_factor * target_weight + 1, - *information_sets.get_overlapping_family_policy( + information_sets.get_overlapping_family_policy( self._overlapping_tuning, len(self._basis) ), ) diff --git a/src/qldpc/codes/distance/_search/information_sets.py b/src/qldpc/codes/distance/_search/information_sets.py index 693d87704..5852ba871 100644 --- a/src/qldpc/codes/distance/_search/information_sets.py +++ b/src/qldpc/codes/distance/_search/information_sets.py @@ -13,11 +13,11 @@ that a basis family has not yet visited. The coset searches use the disjoint information sets, and for costly queries the families of :func:`get_maximum_disjoint_family` and :func:`get_overlapping_family` (see ``coset_routes.py``). The distance, counting, and coset scans -all enumerate a layer with :func:`advance_support`. The distance and coset searches plan their -steps with :func:`plan_steps`: where a search stops, which is the first step whose bound exceeds a -limit, and how much work the steps up to it take. A distance search whose certificate solves linear -programs instead finds where it stops in each layer by checking the steps in order (see -``Certificate.plan``). +all enumerate a layer with :func:`advance_support`. The distance and coset searches find where +they stop with :func:`find_stopping_step`: the first step whose bound exceeds a limit. +:func:`plan_steps` also adds up the work of the steps up to it. A distance search whose certificate +solves linear programs instead finds where it stops in each layer by checking the steps in order +(see ``Certificate.plan``). A :class:`BasisFamily` stacks its information sets into arrays, so that one set of functions serves every caller, in the subset of Python that numba compiles, and some of them are compiled (see @@ -88,9 +88,13 @@ class Certificate: once the first sets in family order have completed it, in whatever order their words were visited. A set that is interrupted within a layer does not advance either bound. The bounds do not bound the weights of witnesses already found, and they do not transfer to another basis or - CSS sector. Each bound is the coverage bound: the occupancy of the sets (see - ``get_step_lower_bounds``) divided by ``max_coverage``, rounded up, and then rounded up to a - multiple of ``weight_divisor``. With ``pivot_blocks``, each bound is the larger of the coverage + CSS sector. A search that skips supports by completion tables completes a layer without + visiting them, and then the bounds hold only for the eligible words lighter than its best weight + (Section 3.9 of ``brouwer_zimmermann.md``). + + Each bound is the coverage bound: the occupancy of the sets (see ``get_step_lower_bounds``) + divided by ``max_coverage``, rounded up, and then rounded up to a multiple of + ``weight_divisor``. With ``pivot_blocks``, each bound is the larger of the coverage bound and the bound of ``get_lp_lower_bound``, which can be larger when coverage is uneven. Attributes: @@ -169,19 +173,15 @@ def plan(self, limit: int, layer: int) -> int: The search stops at its stopping step: the first step of the layer whose ``step_lower_bound`` exceeds ``limit``, or the last step, ``(layer, len(ranks))``, if none does. Returns the number of sets that have completed the layer at that step. Without - linear programs, ``plan_steps`` finds the step. A certificate that solves linear programs - checks the steps in order with ``certifies``, except the last step, which stops the search - anyway. So it solves a program only for a step before the stopping step or at it, never for - the last step, and only when the coverage bound of the step does not exceed the limit. + linear programs, ``find_stopping_step`` finds the step. A certificate that solves linear + programs checks the steps in order with ``certifies``, except the last step, which stops the + search anyway. So it solves a program only for a step before the stopping step or at it, + never for the last step, and only when the coverage bound of the step does not exceed the + limit. """ num_sets = len(self.ranks) if not self._solves_linear_programs: - # Empty binomials ask ``plan_steps`` for the step only, and a plan that leaves out no - # sets reads no ranks, so both are empty placeholders. - no_values = np.empty(0, dtype=np.int64) - _, num_completed, _ = plan_steps( - self._step_lower_bounds, limit, layer, layer, no_values, no_values, False - ) + _, num_completed = find_stopping_step(self._step_lower_bounds, limit, layer, layer) return num_completed for num_completed in range(1, num_sets): if self.certifies(layer, num_completed, limit + 1): @@ -195,8 +195,8 @@ def get_num_sets_to_certify(self, weight: int, upper_bound: int) -> int: unvisited eligible word weighs at least ``upper_bound``. Returns the number of sets if no fewer sets suffice, in which case the bound may need later layers too. The number is that of the stopping step of ``plan`` in the layer for the limit ``upper_bound - 1``. Without - linear programs, ``plan_steps`` finds that step, so the distance and coset searches find - their stopping steps by one shared rule. + linear programs, ``find_stopping_step`` finds that step, so the distance and coset searches + find their stopping steps by one shared rule. """ return self.plan(upper_bound - 1, weight) @@ -379,33 +379,49 @@ def prune_family(basis_family: BasisFamily, upper_bound: int) -> BasisFamily: ) +class OverlappingFamilyPolicy(NamedTuple): + """How stage three builds a replacement basis family (see ``get_overlapping_family``). + + Searches take it from :func:`get_overlapping_family_policy`. + + Attributes: + max_sets: The most sets of a replacement. + minimum_work: The estimated work below which a basis family is not replaced. + attempts_per_set: The column orders that the construction tries for each set that it may + add. + binomials: ``get_capped_binomials(dimension)``, for the dimension of the basis. + rng: The generator that breaks ties between equally used columns. + """ + + max_sets: int + minimum_work: int + attempts_per_set: int + binomials: npt.NDArray[np.int64] + rng: np.random.Generator + + @compiled_on_request def get_overlapping_family( basis: npt.NDArray[np.uint8], labels: npt.NDArray[np.uint8], basis_family: BasisFamily, upper_bound: int, - max_sets: int, - minimum_work: int, - attempts_per_set: int, - binomials: npt.NDArray[np.int64], - rng: np.random.Generator, + policy: OverlappingFamilyPolicy, ) -> BasisFamily: """Stage three: build a replacement basis family of overlapping full sets that does less work. ``basis_family`` holds the full and partial sets of one independent binary basis, whose quotient labels ``labels`` have no columns if it has none. A replacement holds only full sets, with their labels transformed alongside their rows, and is certified with its maximum coverage of a - Hamming coordinate. ``max_sets`` bounds its size and the effort spent building it: more sets - can avoid costly coefficient layers in hard searches, while fewer reduce set-up work. No - replacement is built if the estimated work of ``basis_family`` is below ``minimum_work``. The - construction tries ``attempts_per_set`` column orders for each set that it may add. - - Returns the replacement, which has no sets if no replacement does less work. ``binomials`` are - ``get_capped_binomials(dimension)``, and ``rng`` breaks ties between equally used columns. - Searches take ``max_sets``, ``minimum_work``, ``attempts_per_set``, ``binomials``, and ``rng`` - from :func:`get_overlapping_family_policy`. + Hamming coordinate. ``max_sets`` of ``policy`` bounds its size and the effort spent building + it: more sets can avoid costly coefficient layers in hard searches, while fewer reduce set-up + work. No replacement is built if the estimated work of ``basis_family`` is below + ``minimum_work``. The construction tries ``attempts_per_set`` column orders for each set that + it may add. + + Returns the replacement, which has no sets if no replacement does less work. """ + max_sets, minimum_work, attempts_per_set, binomials, rng = policy dimension = basis.shape[0] columns = np.flatnonzero(basis.sum(axis=0)) num_columns = len(columns) @@ -451,8 +467,9 @@ def get_overlapping_family( best_num_sets = 0 best_coverage = 1 packed_basis, packed_labels = pack_bits(basis), pack_bits(labels) - # Saturate the number of attempts, which compiled code computes with 64-bit integers. - num_open = max(max_sets - num_full, 1) + # Saturate the number of attempts, which compiled code computes with 64-bit integers. At least + # one set is open, since the family has fewer full sets than the limit. + num_open = max_sets - num_full num_attempts = ( attempts_per_set * num_open if attempts_per_set <= _tuning.MAX_VALUE // num_open @@ -510,17 +527,15 @@ def get_overlapping_family( ) -def get_overlapping_family_policy( - tuning: SearchTuning, dimension: int -) -> tuple[int, int, int, npt.NDArray[np.int64], np.random.Generator]: - """The arguments of stage three after ``upper_bound``, from ``max_sets`` on, for a search. +def get_overlapping_family_policy(tuning: SearchTuning, dimension: int) -> OverlappingFamilyPolicy: + """The policy of stage three for a search of a basis of a given dimension. - They are the limit, the minimum work, and the attempts of ``tuning``, the capped binomials of + It holds the limit, the minimum work, and the attempts of ``tuning``, the capped binomials of the dimension, and a new generator with a fixed seed. So every search, of the distance, of the counts, or of cosets (see ``coset_routes.py``), chooses the same replacement for the same basis family and tuning. """ - return ( + return OverlappingFamilyPolicy( tuning.max_information_sets, tuning.overlapping_minimum_work, tuning.overlapping_attempts_per_set, @@ -932,39 +947,37 @@ def plan_steps( ) -> tuple[int, int, int]: """Plan a search of layers ``first_layer`` to ``last_layer`` that stops above a limit. - The coset searches find where they stop with this function, and so does the distance search, - through ``Certificate.plan``, when its certificate solves no linear programs. The work - estimates of the distance search add up the work of its steps with it. The table - ``step_bounds`` has a row for each layer from zero to the dimension, and a column for each set - after column zero. The steps of layer ``w`` are ``(w, 1)`` to ``(w, num_sets)``, in that order - and layer by layer, and the bound of step ``(w, j)`` is ``step_bounds[w, j]``, as in the table - of ``get_step_lower_bounds``. These bounds must not decrease in the order of the steps, which - holds for that table and for any nondecreasing function of it, such as the tables of - ``Certificate.get_weight_bounds``. + The work estimates of the distance search and of the routes of the coset searches add up the + work of the steps of a search with this function, and searches that need only where they stop + call ``find_stopping_step``. The table ``step_bounds`` has a row for each layer from zero to + the dimension, and a column for each set after column zero. The steps of layer ``w`` are + ``(w, 1)`` to ``(w, num_sets)``, in that order and layer by layer, and the bound of step + ``(w, j)`` is ``step_bounds[w, j]``, as in the table of ``get_step_lower_bounds``. These bounds + must not decrease in the order of the steps, which holds for that table and for any + nondecreasing function of it, such as the tables of ``Certificate.get_weight_bounds``. The stopping step ``(w*, j*)`` is the first step of the layers whose bound exceeds ``limit``, - which a binary search finds, or ``(last_layer, num_sets)`` if there is none. At that step, set - ``i`` has completed layer ``w_i = w*`` if ``i < j*``, and ``w* - 1`` otherwise. A set of rank - deficiency ``s`` adds to the bound of the step only if ``w_i >= s`` (see + or ``(last_layer, num_sets)`` if there is none, which ``find_stopping_step`` finds. At that + step, set ``i`` has completed layer ``w_i = w*`` if ``i < j*``, and ``w* - 1`` otherwise. A set + of rank deficiency ``s`` adds to the bound of the step only if ``w_i >= s`` (see ``get_step_lower_bounds``). The other sets are the unused sets of the step. If ``skip_unused`` is True, the plan leaves them out, as a coset search skips them, and reads the ranks of the sets from ``ranks``. Otherwise, it reads no ranks. The work adds ``binomials[w]`` candidate/set evaluations for each set that the plan keeps and - each layer ``w`` from ``first_layer`` to ``w_i``. Empty ``binomials`` ask for the step only, - with zero work. The work saturates at ``_WORK_CAP``: it is the cap once a layer up to ``w*`` - has a binomial above ``_WORK_CAP // num_sets``, or once the sum reaches the cap after a layer. - With the binomials of ``get_capped_binomials``, a plan from layer one for the limit ``U - 1`` - gives the work estimate of a distance search (Section 3.6 of ``brouwer_zimmermann.md``). A plan - from layer zero for a maximum weight ``max_weight`` that skips unused sets gives the number of - words that a coset search weighs, or skips by completion tables (Section 3.9 of - ``brouwer_zimmermann.md``), after the word of its syndrome, if no word of the coset weighs at - most ``max_weight``. + each layer ``w`` from ``first_layer`` to ``w_i``. Empty ``binomials`` give zero work. The work + saturates at ``_WORK_CAP``: it is the cap once a layer up to ``w*`` has a binomial above + ``_WORK_CAP // num_sets``, or once the sum reaches the cap after a layer. With the binomials of + ``get_capped_binomials``, a plan from layer one for the limit ``U - 1`` gives the work estimate + of a distance search (Section 3.6 of ``brouwer_zimmermann.md``). A plan from layer zero for a + maximum weight ``max_weight`` that skips unused sets gives the number of words that a coset + search weighs, or skips by completion tables (Section 3.9 of ``brouwer_zimmermann.md``), after + the word of its syndrome, if no word of the coset weighs at most ``max_weight``. Returns: The stopping step ``(w*, j*)``, and the work. """ - layer, num_completed = _find_stopping_step(step_bounds, limit, first_layer, last_layer) + layer, num_completed = find_stopping_step(step_bounds, limit, first_layer, last_layer) work = _get_steps_work( step_bounds.shape[0] - 1, step_bounds.shape[1] - 1, @@ -978,6 +991,31 @@ def plan_steps( return layer, num_completed, work +@always_compiled +def find_stopping_step( + step_bounds: npt.NDArray[np.int64], limit: int, first_layer: int, last_layer: int +) -> tuple[int, int]: + """Find the stopping step of a search of layers ``first_layer`` to ``last_layer``. + + The steps and their bounds are those of ``plan_steps``, and so is the stopping step: the first + step of the layers whose bound exceeds ``limit``, or ``(last_layer, num_sets)`` if there is + none. The bounds do not decrease in the order of the steps, so a binary search finds it. The + coset searches find where they stop with this function, and so does the distance search, + through ``Certificate.plan``, when its certificate solves no linear programs. + """ + num_sets = step_bounds.shape[1] - 1 + low, high = 0, (last_layer - first_layer + 1) * num_sets + while low < high: + middle = (low + high) // 2 + if step_bounds[first_layer + middle // num_sets, 1 + middle % num_sets] > limit: + high = middle + else: + low = middle + 1 + if low == (last_layer - first_layer + 1) * num_sets: + return last_layer, num_sets + return first_layer + low // num_sets, 1 + low % num_sets + + @always_compiled(inline=True) def get_completed_layer(set_index: int, layer: int, num_completed: int) -> int: """The last layer that set ``set_index`` has completed at step ``(layer, num_completed)``.""" @@ -1077,24 +1115,6 @@ def _estimate_optimistic_work( return optimistic_work -@always_compiled -def _find_stopping_step( - step_bounds: npt.NDArray[np.int64], limit: int, first_layer: int, last_layer: int -) -> tuple[int, int]: - """Find the stopping step of ``plan_steps`` by binary search.""" - num_sets = step_bounds.shape[1] - 1 - low, high = 0, (last_layer - first_layer + 1) * num_sets - while low < high: - middle = (low + high) // 2 - if step_bounds[first_layer + middle // num_sets, 1 + middle % num_sets] > limit: - high = middle - else: - low = middle + 1 - if low == (last_layer - first_layer + 1) * num_sets: - return last_layer, num_sets - return first_layer + low // num_sets, 1 + low % num_sets - - @always_compiled def _get_steps_work( dimension: int, diff --git a/src/qldpc/codes/distance/_search/information_sets_test.py b/src/qldpc/codes/distance/_search/information_sets_test.py index ac0366c03..2a6013ff9 100644 --- a/src/qldpc/codes/distance/_search/information_sets_test.py +++ b/src/qldpc/codes/distance/_search/information_sets_test.py @@ -714,7 +714,7 @@ def _get_overlapping_family( np.empty((len(basis), 0), dtype=np.uint8) if labels is None else labels, basis_family, upper_bound, - *information_sets.get_overlapping_family_policy(tuning, len(basis)), + information_sets.get_overlapping_family_policy(tuning, len(basis)), ) diff --git a/src/qldpc/codes/distance/_tuning.py b/src/qldpc/codes/distance/_tuning.py index c8c3b98fe..2bfc221f9 100644 --- a/src/qldpc/codes/distance/_tuning.py +++ b/src/qldpc/codes/distance/_tuning.py @@ -21,9 +21,9 @@ - Tuning changes only how long a search takes and how much memory it uses. It never changes an exact distance, a distance limited by ``max_weight``, or a count of logical operators. -- With an early exit (``cutoff`` of two or more), a distance search returns the weight of the first - nontrivial word that it finds at or below ``cutoff``, and which word comes first depends on - tuning. The result is still an upper bound between the distance and ``cutoff``, as documented by +- With an early exit (``cutoff`` of two or more), a distance search returns the weight of a + nontrivial word at or below ``cutoff`` that it finds, and which word that is depends on tuning. + The result is still an upper bound between the distance and ``cutoff``, as documented by :func:`~qldpc.codes.distance.get_distance_quantum`. - Under every tuning, a query of ``CosetLeaderSearch`` decides the same way whether a coset has a word of weight at most the query's maximum weight, and with the connected-support method it gives diff --git a/src/qldpc/codes/distance/_tuning_test.py b/src/qldpc/codes/distance/_tuning_test.py index f9ba72673..4b5ade0bd 100644 --- a/src/qldpc/codes/distance/_tuning_test.py +++ b/src/qldpc/codes/distance/_tuning_test.py @@ -171,6 +171,7 @@ def test_tunings_change_routes_but_not_results(costly_stages: CostlyStages) -> N distances = (fresh.get_distance_exact(Pauli.X), fresh.get_distance_exact(Pauli.Z)) counter_code = codes.CSSCode(code.matrix_x, code.matrix_z) counts = lc.CSSLogicalCounter(counter_code).count_logical_ops(count_pauli=True) + assert counts.weight is not codes.NO_DISTANCE macwilliams = lc.CSSLogicalCounter(counter_code).count_logical_ops( counts.weight, method="macwilliams" ) @@ -282,7 +283,7 @@ def recording(*args: Any, **kwargs: Any) -> Any: record(bz, "enumeration_is_cheaper_than_setup", "max_dimension", "setup_cost_factor") record(information_sets, "_get_disjoint_sets", "random_column_orders") - record(information_sets, "get_overlapping_family", "attempts_per_set") + record(information_sets, "get_overlapping_family_policy", "tuning", prefix="overlapping ") record(information_sets, "get_partitioned_family", "tuning") record(information_sets.BasisFamily, "get_search_certificate", "lp_max_denominator") record(bz, "get_nested_distance_brute_force", "batch_exponent", prefix="bz ") @@ -337,7 +338,7 @@ def prepare(*args: Any, **kwargs: Any) -> Any: "max_dimension": {3}, "setup_cost_factor": {1}, "random_column_orders": {5}, - "attempts_per_set": {3}, + "overlapping tuning": {tuning}, "tuning": {tuning}, "lp_max_denominator": {97}, "bz batch_exponent": {6}, diff --git a/src/qldpc/codes/distance/brouwer_zimmermann.md b/src/qldpc/codes/distance/brouwer_zimmermann.md index 26d2a9e99..b57a08de1 100644 --- a/src/qldpc/codes/distance/brouwer_zimmermann.md +++ b/src/qldpc/codes/distance/brouwer_zimmermann.md @@ -490,7 +490,8 @@ $$ W = m\sum_{w=1}^{t_* - 1} \binom{r}{w} + j_* \binom{r}{t_*}. $$ -qLDPC computes this estimate with `plan_steps`, the function that also finds where a search stops ([Section 4.1](#41-module-map-and-layering)). +qLDPC computes this estimate with `plan_steps`, which finds where the search stops with `find_stopping_step`. +Searches whose certificates solve no linear programs find their stopping steps the same way ([Section 4.1](#41-module-map-and-layering)). The estimate is pessimistic, because it assumes that no lighter witness is found before the certificate reaches `U`. The sum starts at `w = 1`, so it includes the systematic rows. When the comparison with exhaustive enumeration is made, layer one has already been checked, so `W` slightly overestimates the remaining BZ work. @@ -499,7 +500,7 @@ The estimate leaves out the search bases that the query skips. If no word of the coset weighs at most `ω`, the query examines, or skips by completion tables ([Section 3.9](#39-skipping-completions-by-completion-tables)), exactly that many words besides the word of its syndrome, and otherwise at most that many. Actual runtime also depends on basis construction, the width of packed words, and early improvements of the witness, so `W` is not a runtime guarantee. -Exhaustive enumeration is chosen at three points, which [Section 4.2](#42-flow-and-exit-points) lists. +Exhaustive enumeration is chosen at two points, which [Section 4.2](#42-flow-and-exit-points) lists. It enumerates words in vectorized, packed blocks and updates them with Gray-code steps, and it does not use the BZ certificate. ### 3.7 Images of an information set under automorphisms @@ -602,14 +603,15 @@ L(4)=\left\lceil\frac{3(4+1)}{2}\right\rceil=8, $$ which reaches eight even before rounding to an even weight. -The search stops after layer four. +qLDPC stops even sooner, part way through layer four ([Section 2.5](#25-the-stopping-certificate)): after the first search basis of that layer, the bound is $\lceil(5+4+4)/2\rceil = 7$, which rounds up to eight. +The search therefore scans one search basis in layer four. | Family | Last layer | Candidate/set evaluations in layers two and above | | --- | ---: | ---: | | Disjoint, one search basis | 6 | $\sum_{w=2}^{6}\binom{41}{w} = 5{,}358{,}536$ | -| Overlapping, three search bases, `C = 2` | 4 | $3\sum_{w=2}^{4}\binom{41}{w} = 338{,}250$ | +| Overlapping, three search bases, `C = 2` | 4 | $3\sum_{w=2}^{3}\binom{41}{w} + \binom{41}{4} = 135{,}710$ | -The overlapping family needs fewer evaluations because it avoids layers five and six. +The overlapping family needs fewer evaluations because it avoids layers five and six, and most of layer four. The counts cover only enumeration; building and packing the search bases costs extra. ### 3.9 Skipping completions by completion tables @@ -715,7 +717,8 @@ A search for the minimum weight of a coset skips at `h(ω)` too until it finds a #### When a search uses the tables -Every layer of two or more rows of a distance search or count that scans one information set uses `scan_information_set` of [_search/affine_search.py](_search/affine_search.py), a scan specialized to that case, with or without tables. +Every layer of two or more rows of a distance search that scans one information set uses `scan_information_set` of [_search/affine_search.py](_search/affine_search.py), a scan specialized to that case, with or without tables. +A count scans with its own kernels in [_counting/scan.py](_counting/scan.py), which use the tables as described below. Tables of a basis with `r = 150` rows and 138 tail columns, in blocks of sixteen columns, take 75 MiB, and about a twentieth of a second to build on the machine described below. That is worthwhile for a layer of `10^10` supports, but not for an easy search. A search considers tables for a layer of at least `completion_bound_minimum_layer_work` supports ([Section 4.4](#44-heuristics-and-cost-models)). @@ -831,8 +834,8 @@ Distance searches and counts use [_search/completion_bounds.py](_search/completi The distance search and `CosetLeaderSearch` scan their layers with the compiled entry points of [_search/affine_search.py](_search/affine_search.py), which run one body, `_scan`, with constant flags. `scan_family` scans every set of a family in lockstep, `scan_information_set` scans one information set with or without completion tables, and `scan_coset_layer` and `scan_coset_layer_by_css_weight` scan one layer of one set of a coset and keep its lightest word. `scan_coset_layer_by_tails` scans one layer of one set of a coset by the tails of its words, without the pivot columns of the set ([Section 3.5](#35-packed-enumeration-with-pivot-column-deletion)), and keeps the support of its lightest word. -`scan_coset_layer_by_tails_with_tables` does the same for a layer of two or more rows, and skips the completions of a prefix by tables if their bound reaches a skip bound ([Section 3.9](#39-skipping-completions-by-completion-tables)). -`scan_coset_layer_by_tails_with_tables_below_best` skips them instead if their bound reaches the best weight so far, for the searches for minimum weights of [Section 3.10](#310-minimum-weights-of-cosets). +`scan_coset_layer_by_tails_with_tables` does the same for a layer of two or more rows, and skips the completions of a prefix by tables if their bound reaches the lesser of a skip bound and the best weight so far ([Section 3.9](#39-skipping-completions-by-completion-tables)). +So a query skips at its skip bound, since its scans return before their best weight falls below it, and a search for minimum weights of [Section 3.10](#310-minimum-weights-of-cosets) skips at the best weight, which never exceeds its skip bound. Numba compiles the body separately for each entry point, so no flag costs a branch for each word. Each search builds an `AffineInstance` of [_search/affine_instance.py](_search/affine_instance.py) for its problem before it scans, rather than for each layer or query. The instance holds the packed rows, labels, and pivots of the sets that the search scans, the ranks of those sets, and the certificate of their steps. @@ -857,8 +860,9 @@ A query of a coset search asks the same of a coset of its code, whose words have It starts from the **word of its syndrome**, a word of the coset that sums the rows of a right inverse of the checks for the set bits of the syndrome. It stops at its **stopping step**, the first step whose bound exceeds the maximum weight, unless it finds a light enough word sooner. It skips each rank-deficient search basis whose deficiency keeps it out of the bound of that step. -`plan_steps` of [_search/information_sets.py](_search/information_sets.py) plans the steps of the coset searches and of the distance search. -It finds the stopping step by a binary search of steps whose bounds never decrease, and it adds up the work of the steps up to the stopping step ([Section 3.6](#36-exhaustive-enumeration-as-a-fallback)). +`find_stopping_step` of [_search/information_sets.py](_search/information_sets.py) finds the stopping steps of the coset searches, and of the distance searches whose certificates solve no linear programs, by a binary search of steps whose bounds never decrease. +A certificate that solves linear programs checks steps in order instead, through `Certificate.plan` (see below). +`plan_steps` also adds up the work of the steps up to the stopping step ([Section 3.6](#36-exhaustive-enumeration-as-a-fallback)). Each family of a coset search keeps the bounds of its steps in a **planning buffer**, with a table for each parity of the cosets if the search rounds to parity, and one table otherwise. A query takes the table of the parity of the word of its syndrome. Each table has one more row that caches a stopping step. @@ -869,24 +873,25 @@ Otherwise, a query takes the tables of the **route** of its maximum weight, whic A route is a family of search bases of the code, with its coset instance and planning buffer: the disjoint sets, or a family that the routes build for costly queries. `find_first_not_within_max_weight` answers a batch with the route of its maximum weight, which it takes once, and the batch of the X parts and that of the Z parts each take their own. `_get_min_weights` of `CosetLeaderSearch` finds the minimum weights of [Section 3.10](#310-minimum-weights-of-cosets) for the residual-weight estimator. -`_CosetSearch.get_min_weights` answers a batch on the route of the maximum weight ([Section 4.4](#44-heuristics-and-cost-models)), with `_find_min_weights`, or with `_find_min_weights_by_tails` on a route by tails. -These searches find their stopping steps with `plan_steps`, and leave the cached steps of the planning buffer to the queries. +`_CosetSearch.get_min_weights` answers a batch on the route of the maximum weight ([Section 4.4](#44-heuristics-and-cost-models)), with `_find_min_weights`. +These searches find their stopping steps with `find_stopping_step`, and leave the cached steps of the planning buffer to the queries. The gate tests the predicted examinations of a query with the disjoint family against its threshold, to decide whether the query may take another family ([Section 4.4](#44-heuristics-and-cost-models)). A route is a **route by tails** if those with the route's own family, which can do less work, also reach the threshold. -The queries of a route by tails weigh words by their tails ([Section 3.5](#35-packed-enumeration-with-pivot-column-deletion)), and run `_search_coset_by_tails` instead of `_search_coset`. +The queries of a route by tails weigh words by their tails ([Section 3.5](#35-packed-enumeration-with-pivot-column-deletion)): they pass the record of tails of the route to the compiled searches, such as `_search_coset`, in place of None. +Numba compiles each search separately for None and for a record of tails, so the choice costs nothing for each word. Every route by tails of a family gives the latest record `_CosetTails` of the family. The record holds the tails of the rows of each search basis, and of the rows of the right inverse of the checks, from which a query adds up the tail of each offset. It also holds the completion tables of the search bases that the queries of the route scan in layers of two or more rows, if they fit in `completion_bound_max_table_bytes`. The routes build the tails of a family once, and the tables of each of its search bases once, for the first route that needs them (`_get_coset_tails`, `completion_bounds.considers_tables`, and `completion_bounds.build_tables_within`). They keep only the latest record of a family, so they hold the tables of each search basis once. -`_search_coset_by_tails` scans each layer of two or more rows of a search basis that has tables with `scan_coset_layer_by_tails_with_tables`. +A query by tails scans each layer of two or more rows of a search basis that has tables with `scan_coset_layer_by_tails_with_tables`. It scans the other layers with `scan_coset_layer_by_tails`. -It visits the words of `_search_coset` in the same order, except those that it skips, so it finds the same first word within the maximum weight. +It visits the words of a query of whole words in the same order, except those that it skips, so it finds the same first word within the maximum weight. Without one, it can return another word and another lower bound ([Section 3.9](#39-skipping-completions-by-completion-tables)). It keeps the search basis and support of its lightest word, from which it rebuilds the word once, with the full rows. The search of Pauli errors by CSS weight has no routes by tails, since weighing a word by its CSS weight needs to know which deleted pivot columns lie in the X part. The distance search finds the stopping step of each layer with `Certificate.plan`, through `Certificate.get_num_sets_to_certify`. -`Certificate.plan` calls `plan_steps` when its certificate solves no linear programs. +`Certificate.plan` calls `find_stopping_step` when its certificate solves no linear programs. A certificate with linear programs instead checks the steps of the layer but the last in order with `certifies`, which solves each program only when the coverage bound falls short. For the CSS weight `max(|x|, |z|)` on a CSS state, whose stabilizer group is generated by stabilizers with no Z part and stabilizers with no X part, `CosetLeaderSearch` instead runs this search on two classical codes, spanned by the X parts of the first kind and by the Z parts of the second, with words weighed by their Hamming weight. @@ -902,7 +907,7 @@ An integer returned by preparation can incorporate a witness from another sector | Public route in [core.py](core.py) | Outcome | | --- | --- | -| Trivial input | Returns `n+1`, capped at `max_weight+1`, for an empty classical code, a quantum problem with no logical rows, or `V = S` | +| Trivial input | Returns `NO_DISTANCE`, whatever `max_weight` is, for an empty classical code, a quantum problem with no logical rows, or `V = S` | | `method="brute_force"` | Runs exhaustive enumeration over an independent basis, adapted to `S ⊆ V` for a quantum problem | `_prepare_brouwer_zimmermann_search` considers the following exits in order. @@ -919,7 +924,7 @@ The disjoint family, its pruning, and any overlapping or partitioned replacement | Exhaustive by work | Runs exhaustive enumeration if the tuning field `allow_enumeration_by_work` is true and `_enumeration_is_cheaper_than_search` finds that the eligible-word count is at most the estimated BZ work of [Section 3.6](#36-exhaustive-enumeration-as-a-fallback) | | Ready search | Returns a `_BrouwerZimmermannSearch` that scans from layer two | -Exhaustive enumeration is chosen automatically at three points: the direct shortcut, preparation's size check, and preparation's work check. +Exhaustive enumeration is chosen automatically at two points: preparation's size check, and preparation's work check. | Search location | Outcome | | --- | --- | @@ -927,7 +932,7 @@ Exhaustive enumeration is chosen automatically at three points: the direct short | `_BrouwerZimmermannSearch.finished` | Marks a search finished if `U ≤ cutoff`, `L ≥ U`, or every layer has been enumerated (`next_weight > dimension`) | | `get_distance_brouwer_zimmermann_many` | Returns the shared best bound when it is at most `cutoff` after an advance, or when no search can advance | -The tests `test_brouwer_zimmermann_setup_exits`, `test_brouwer_zimmermann_setup_cutoffs`, `test_rows_of_an_automorphism_family_can_reach_the_cutoff`, and `test_brouwer_zimmermann_exhaustive_dispatch` exercise these exits. +The tests `test_brouwer_zimmermann_setup_exits`, `test_brouwer_zimmermann_setup_cutoffs`, `test_rows_of_an_automorphism_family_can_reach_the_cutoff`, and `test_enumeration_by_work_can_be_forbidden` exercise these exits. ### 4.3 Concept-to-code glossary @@ -935,7 +940,7 @@ The tests `test_brouwer_zimmermann_setup_exits`, `test_brouwer_zimmermann_setup_ | --- | --- | | Upper bound `U` | `best`, or `upper_bound` as a parameter | | Lower certificate `L` | `lower_bound` of a search, computed by `Certificate.lower_bound(completed_weight)` from the certificate's `dimension`, `ranks`, `weight_divisor`, and `max_coverage`; it bounds the eligible words lighter than `U`, so a bound on the distance is `min(L, U)` | -| Certificate within a layer `L(t, j)` | `Certificate.step_lower_bound(weight, num_completed)` and the table of `get_step_lower_bounds`. `Certificate.get_num_sets_to_certify` gives the smallest `j` that certifies a bound. `plan_steps` gives the first step whose bound exceeds a limit, the stopping step, and the work of the steps up to it. `Certificate.plan` gives the stopping step of one layer, with linear programs if the certificate solves them. `Certificate.get_weight_bounds` divides the table by a weight factor, rounding up, after rounding it up to a parity if it is given one. A coset search keeps the results in its planning buffer, `planning` of `_CosetSearchTables`: a table for each parity of its cosets if it rounds to parity, and one table otherwise, each with a cached stopping step. | +| Certificate within a layer `L(t, j)` | `Certificate.step_lower_bound(weight, num_completed)` and the table of `get_step_lower_bounds`. `Certificate.get_num_sets_to_certify` gives the smallest `j` that certifies a bound. `find_stopping_step` gives the first step whose bound exceeds a limit, the stopping step, and `plan_steps` also the work of the steps up to it. `Certificate.plan` gives the stopping step of one layer, with linear programs if the certificate solves them. `Certificate.get_weight_bounds` divides the table by a weight factor, rounding up, after rounding it up to a parity if it is given one. A coset search keeps the results in its planning buffer, `planning` of `_CosetSearchTables`: a table for each parity of its cosets if it rounds to parity, and one table otherwise, each with a cached stopping step. | | Certificate from linear programming | `get_lp_lower_bound`, which a `Certificate` with `pivot_blocks` combines with the coverage-divided certificate; `BasisFamily.get_search_certificate` adds the blocks to costly searches of overlapping families | | Family of search bases | `BasisFamily`, which stacks the search bases into arrays and gives their `Certificate`. A distance search keeps the `AffineInstance` of [_search/affine_instance.py](_search/affine_instance.py) that packs them, with the certificate. A coset search keeps the generators, pivots, and planning buffer of each family of its routes in `_CosetSearchTables`, and the tails and completion tables of the family of each route by tails in `_CosetTails`. | | Automorphisms, generators of `G` | `automorphisms`, whose row `i` maps column `j` to column `automorphisms[i, j]`, as returned by `get_automorphisms`; `_found_qudit_automorphisms` and `_get_qudit_automorphisms` of `QuditCode` | @@ -977,8 +982,8 @@ The Tuning section of [README.md](README.md) explains how searches read the thre A search basis whose block has rank deficiency `d = r - r_i` has a term `b_i(t)` that is zero for `t < d`. The function removes such a search basis when the blocks that are information sets reach `U` by layer `d`, so the removal leaves the result unchanged and saves enumeration. It removes nothing if no block is an information set. -- **`get_overlapping_family`** tries to lower the work estimate of the disjoint family, with the thresholds `max_information_sets` and `overlapping_minimum_work`, which it receives as `max_sets` and `minimum_work`. - It returns a basis family with no sets if the existing work is below `minimum_work`, if the blocks already use every nonzero column, if the family is at `max_sets`, or if even the optimistic coverage `ceil(m*r/N)` cannot reduce the estimate. +- **`get_overlapping_family`** tries to lower the work estimate of the disjoint family, with the thresholds `max_information_sets` and `overlapping_minimum_work`, which it receives as the fields `max_sets` and `minimum_work` of its `OverlappingFamilyPolicy`. + It returns a basis family with no sets if no block is an information set, if the existing work is below `minimum_work`, if the blocks already use every nonzero column, if the family is at `max_sets`, or if even the optimistic coverage `ceil(m*r/N)` cannot reduce the estimate. Otherwise it builds search bases by row reducing on column orders that put the least-used columns first, breaking ties with a generator that `get_overlapping_family_policy` seeds afresh for every call, and keeps the cheapest family built. It makes `overlapping_attempts_per_set` attempts for each set that it may add, and a repeated pivot set is discarded but still consumes an attempt. - **`get_partitioned_family`** builds families of information sets in which no column appears more than `b` times, for `b` from two to `max_partition_coverage`. @@ -1015,7 +1020,8 @@ The Tuning section of [README.md](README.md) explains how searches read the thre Its cost model was fitted to the scan of a distance search; counting weighs full rows and labels, so its scans cost more for each support, and the model then predicts conservatively. - **`Certificate.work`, `estimate_automorphism_family_work`, `_enumeration_is_cheaper_than_search`, and `enumeration_is_cheaper_than_setup`** are cost models, which Sections [3.6](#36-exhaustive-enumeration-as-a-fallback) and [3.7](#37-images-of-an-information-set-under-automorphisms) explain. They are not runtime guarantees. - `Certificate.work` and `estimate_automorphism_family_work` are the work of plans of `plan_steps`, the function that also finds where searches stop ([Section 4.1](#41-module-map-and-layering)). + `Certificate.work` and `estimate_automorphism_family_work` are the work of plans of `plan_steps`, which finds where searches stop with `find_stopping_step`. + Searches whose certificates solve no linear programs find their stopping steps the same way ([Section 4.1](#41-module-map-and-layering)). `enumeration_is_cheaper_than_setup` compares the setup cost of BZ with packed enumeration of all words, with the fields `enumeration_max_dimension` and `enumeration_setup_cost_factor`. ### 4.5 Behavior spanning files @@ -1034,19 +1040,19 @@ The Tuning section of [README.md](README.md) explains how searches read the thre The limit also caps `cutoff` at `max_weight`. Without that cap, a cutoff at or above `max_weight+1` would end the search at once on the assumed witness, even if a lighter word exists. Passing `cutoff=max_weight` as well stops the search at the first word within the limit. - Inputs without nontrivial words have no distance and give `n+1`, capped at `max_weight+1`. - No word on `n` bits or qubits has weight `n+1`. + Inputs without nontrivial words have no distance and give `NO_DISTANCE`, a float `nan` (see `NoDistance` in `core.py`), whatever `max_weight` is. + Inside the searches, a problem without nontrivial words counts as `n+1`, which no word on `n` bits or qubits weighs, so that it never lowers the minimum over problems. For non-homogeneous input, the BZ route doubles the bound along with `cutoff`. - **Units.** Homogeneous input represents a pure-X or pure-Z problem with one support bit per qubit. Non-homogeneous input (`homogeneous=False`) uses general symplectic Pauli rows, as in [Section 2.7](#27-pauli-weight-as-hamming-weight). For non-homogeneous input, the BZ route maps rows to Hamming coordinates, doubles `cutoff`, passes `weight_divisor=2`, and halves the result. The exhaustive routes inside preparation receive the same Hamming-coordinate basis, so they use the doubled `cutoff`. - Only the direct shortcut and `method="brute_force"` use native symplectic weights and the original `cutoff`. + Only `method="brute_force"` uses native symplectic weights and the original `cutoff`. - **Independence.** Callers need not supply independent generators. The classical path uses `get_independent_rows`, and the quantum path uses `get_nested_code_basis`. - The direct shortcut verifies independence itself, so the packed enumerators always receive independent bases. + So the packed enumerators always receive independent bases. - **Resumable layers.** The compiled scan returns after at most `batch_size` supports. A scan of one information set with completion tables ([Section 3.9](#39-skipping-completions-by-completion-tables)) returns after at most `batch_size` supports and lookups together, so it stays responsive while it skips most supports. diff --git a/src/qldpc/codes/distance/core.py b/src/qldpc/codes/distance/core.py index 3febdc814..f4526fc50 100644 --- a/src/qldpc/codes/distance/core.py +++ b/src/qldpc/codes/distance/core.py @@ -5,7 +5,9 @@ from __future__ import annotations import dataclasses +import enum import functools +import math import operator from collections.abc import Callable, Sequence from typing import Literal, SupportsIndex, get_args @@ -37,6 +39,37 @@ "brouwer_zimmermann", "brute_force", "information_sets", "connected_support" ] + +class NoDistance(float, enum.ReprEnum): + """The distance of an object with no nontrivial logical operator to minimize over. + + A code that encodes nothing, or a circuit with no fault that flips a logical observable without + flipping a detector, has no distance. qLDPC then returns ``NO_DISTANCE``, the only member of + this enum, which behaves as a float ``nan``: + + - It prints as ``nan``, also in the representation of a tuple. + - It equals nothing, not even itself, and ordering comparisons with it are False. + - Arithmetic with it gives ``nan``, and NumPy and pandas treat it as missing. + + Test for it with ``distance is NO_DISTANCE``, or with ``math.isnan(distance)``, which also + catches ``nan`` from arithmetic. A result that can be ``NO_DISTANCE`` is annotated + ``int | Literal[NoDistance.NAN]``, so that type checkers treat it as an ``int`` once + ``distance is NO_DISTANCE`` is ruled out. + """ + + NAN = math.nan + + def __repr__(self) -> str: + return float.__repr__(self) + + @classmethod + def _missing_(cls, value: object) -> NoDistance | None: + # Pickling and copying look the member up by its value, which differs from itself. + return cls.NAN if isinstance(value, float) and math.isnan(value) else None + + +NO_DISTANCE = NoDistance.NAN + _CanonicalDistanceMethod = Literal["brouwer_zimmermann", "brute_force", "connected_support"] # The canonical name of every accepted exact-distance method name, including aliases. @@ -59,12 +92,12 @@ def get_distance_classical( block_size: int | None = None, # DEPRECATED: see qldpc._deprecations method: DistanceMethod = "brouwer_zimmermann", max_weight: SupportsIndex | None = None, -) -> int: +) -> int | Literal[NoDistance.NAN]: """Distance of a classical linear binary code. Args: generators: The generator matrix of the classical code whose distance we want to compute. - Both methods reduce redundant rows to a linearly independent basis. + Both available methods reduce redundant rows to a linearly independent basis. cutoff: Exit early and return once an upper bound on distance falls to or below this cutoff. block_size: Deprecated. Use ``qldpc.codes.distance._tuning.use_tuning`` to set ``enumeration_batch_exponent`` and ``information_set_batch_exponent`` instead. @@ -78,16 +111,18 @@ def get_distance_classical( then the smaller of the distance and ``max_weight + 1``. The search stops as soon as it proves that every nontrivial code word is heavier than ``max_weight``, which can be much faster than computing a large distance exactly. The limit also caps ``cutoff``; pass - ``cutoff=max_weight`` as well to stop at the first code word within the limit. A code - without nontrivial code words gives ``n + 1``, capped likewise. + ``cutoff=max_weight`` as well to stop at the first code word within the limit. Returns: The minimum Hamming distance between different code words, or equivalently the minimum - Hamming weight of a nontrivial code word. A code of length ``n`` without nontrivial code - words has no distance and gives ``n + 1``. No code word has this weight, so it is not an - observed weight, even if it is at most ``cutoff``. If ``max_weight`` is given, the result - is the smaller of this value and ``max_weight + 1``, so a result of ``max_weight + 1`` means - that every nontrivial code word is heavier than ``max_weight``. + Hamming weight of a nontrivial code word. As soon as the lightest code word seen so far + has weight ``<= cutoff``, the search stops early and returns that weight, an upper bound on + the distance in ``[distance, cutoff]``; pass ``cutoff=0`` to disable the early exit. If + ``max_weight`` is given, the result is the + smaller of this distance and ``max_weight + 1``, so a result of ``max_weight + 1`` means + that every nontrivial code word is heavier than ``max_weight``. A code without nontrivial + code words has no distance and gives ``NO_DISTANCE`` (see ``NoDistance``), whatever + ``cutoff`` and ``max_weight`` are. Raises: TypeError: If ``max_weight`` is not an integer, or is a bool. @@ -99,10 +134,10 @@ def get_distance_classical( _reject_connected_support(method) cutoff, upper_bound = _get_weight_limit(cutoff, max_weight) assert_binary(generators, "generators") - matrix = as_binary_matrix(generators) + matrix = as_binary_matrix(generators, "generators") basis = get_independent_rows(matrix) if not len(basis): - return _apply_upper_bound(matrix.shape[1] + 1, upper_bound) + return NO_DISTANCE if method == "brouwer_zimmermann": # Labels with no columns exclude no subcode. @@ -134,12 +169,13 @@ def get_distance_quantum( homogeneous: bool = False, method: DistanceMethod = "brouwer_zimmermann", max_weight: SupportsIndex | None = None, -) -> int: +) -> int | Literal[NoDistance.NAN]: """Distance of a binary quantum code. Args: logical_ops: Rows that generate the logical space together with the stabilizers. - Both methods remove redundant rows and exclude the stabilizer row space from witnesses. + Both available methods remove redundant rows and exclude the stabilizer row space from + witnesses. stabilizers: A matrix whose rows represent stabilizers of the code. cutoff: Exit early and return once an upper bound on distance falls to or below this cutoff. block_size: Deprecated. Use ``qldpc.codes.distance._tuning.use_tuning`` to set @@ -158,8 +194,7 @@ def get_distance_quantum( proves that every nontrivial logical operator is heavier than ``max_weight``, which can be much faster than computing a large distance exactly. The limit also caps ``cutoff``, so an early exit still returns the weight of an operator within the limit; pass - ``cutoff=max_weight`` as well to stop at the first such operator. Inputs without - nontrivial logical operators give ``n + 1``, capped likewise. + ``cutoff=max_weight`` as well to stop at the first such operator. Returns: The exact minimum weight of a nontrivial logical operator in ``logical_ops`` modulo @@ -168,11 +203,11 @@ def get_distance_quantum( bound on the true distance, in ``[distance, cutoff]``. With the default ``cutoff=1`` this still returns the exact distance for every valid input (an early return can exceed the distance only when ``cutoff >= 2``). Pass ``cutoff=0`` to disable the early exit and force - the exact minimum. Inputs on ``n`` qubits without nontrivial logical operators have no - distance and give ``n + 1``. No operator has this weight, so it is not an observed weight, - even if it is at most ``cutoff``. If ``max_weight`` is given, the result is the smaller of - this value and ``max_weight + 1``, so a result of ``max_weight + 1`` means that every - nontrivial logical operator is heavier than ``max_weight``. + the exact minimum. If ``max_weight`` is given, the result is the smaller of this value and + ``max_weight + 1``, so a result of ``max_weight + 1`` means that every nontrivial logical + operator is heavier than ``max_weight``. Inputs without nontrivial logical operators have + no distance and give ``NO_DISTANCE`` (see ``NoDistance``), whatever ``cutoff`` and + ``max_weight`` are. Raises: TypeError: If ``max_weight`` is not an integer, or is a bool. @@ -213,7 +248,7 @@ def _get_distance_quantum( max_weight: SupportsIndex | None = None, get_automorphisms: Callable[[], npt.NDArray[np.int_] | None] | None = None, tuning: SearchTuning | None = None, -) -> int: +) -> int | Literal[NoDistance.NAN]: """Compute the distance of ``get_distance_quantum``, letting BZ use qubit automorphisms. ``get_automorphisms``, if given, is a function that returns an array of generators of a group @@ -236,16 +271,13 @@ def _get_distance_quantum( physical_size = num_bits if homogeneous else num_bits // 2 if not homogeneous and num_bits % 2: raise ValueError("Symplectic operators must have an even number of columns") - # Without a nontrivial logical operator there is no distance, and the result is n + 1, which - # exceeds the weight of every operator. - no_distance = _apply_upper_bound(physical_size + 1, upper_bound) if not len(logical_matrix): - return no_distance + return NO_DISTANCE # Separate an independent excluded basis from a complement representing the quotient. basis, labels = get_nested_code_basis(logical_matrix, stabilizer_matrix) if labels.shape[1] == 0: - return no_distance + return NO_DISTANCE if method == "brute_force": distance = get_nested_distance_brute_force( basis, @@ -299,9 +331,8 @@ def _get_distance_css_brouwer_zimmermann( Returns: The minimum exact distance across the supplied sectors and the optional upper bound. A - sector without nontrivial logical operators counts as one more than its width, as in - get_distance_quantum, and no operator has that weight. Any other result at or below - ``cutoff`` is an observed upper bound. + sector without nontrivial logical operators counts as one more than its width, which no + operator has. Any other result at or below ``cutoff`` is an observed upper bound. Raises: ValueError: If there are neither sectors nor an upper bound, if the operators of a sector @@ -321,7 +352,7 @@ def _get_distance_css_brouwer_zimmermann( continue # A sector without nontrivial logical operators has no distance. It contributes one more - # than its width, as get_distance_quantum returns, and no operator has that weight. + # than its width, which no operator has. sector_bound = logical_matrix.shape[1] + 1 no_distance_bound = ( sector_bound if no_distance_bound is None else min(no_distance_bound, sector_bound) @@ -377,8 +408,10 @@ def _as_operator_matrices( """ assert_binary(logical_ops, "logical_ops") assert_binary(stabilizers, "stabilizers") - logical_matrix = as_binary_matrix(logical_ops) - return logical_matrix, as_binary_matrix(stabilizers, num_cols=logical_matrix.shape[1]) + logical_matrix = as_binary_matrix(logical_ops, "logical_ops") + return logical_matrix, as_binary_matrix( + stabilizers, "stabilizers", num_cols=logical_matrix.shape[1] + ) def _get_weight_limit(cutoff: int, max_weight: SupportsIndex | None) -> tuple[int, int | None]: diff --git a/src/qldpc/codes/distance/core_test.py b/src/qldpc/codes/distance/core_test.py index 089e81567..2a837c792 100644 --- a/src/qldpc/codes/distance/core_test.py +++ b/src/qldpc/codes/distance/core_test.py @@ -4,7 +4,10 @@ from __future__ import annotations +import copy import functools +import math +import pickle from collections.abc import Callable from typing import cast, get_args from unittest import mock @@ -48,7 +51,7 @@ def test_deprecated_block_size_argument() -> None: expected = qldpc.codes.get_distance_quantum( logical_ops, stabilizers, homogeneous=True, method="brute_force" ) - searches: list[Callable[[], int]] = [ + searches: list[Callable[[], float]] = [ lambda: qldpc.codes.get_distance_classical(generator, block_size=2), lambda: qldpc.codes.get_distance_classical(generator, block_size=2, method="brute_force"), lambda: qldpc.codes.get_distance_quantum( @@ -157,27 +160,31 @@ def test_connected_support_needs_parity_checks() -> None: def test_brouwer_zimmermann_special_inputs() -> None: """Brouwer-Zimmermann validates matrix shapes and preserves quotient edge cases. - Inputs on n bits or qubits without nontrivial words have no distance and give n + 1, whether - they have no logical rows or only logical rows in the stabilizer span. + Inputs without nontrivial words have no distance and give NO_DISTANCE, whether they have no + logical rows or only logical rows in the stabilizer span. """ - assert qldpc.codes.get_distance_classical(np.empty((0, 4), dtype=np.uint8), cutoff=0) == 5 + no_distance = qldpc.codes.NO_DISTANCE + assert qldpc.codes.get_distance_classical(np.empty((0, 4), dtype=np.uint8), cutoff=0) is ( + no_distance + ) assert ( qldpc.codes.get_distance_quantum( np.empty((0, 4), dtype=np.uint8), np.empty((0, 4), dtype=np.uint8), homogeneous=True, ) - == 5 + is no_distance ) for method in ("brouwer_zimmermann", "brute_force"): distance = qldpc.codes.get_distance_quantum( [[1, 0]], [[1, 0]], homogeneous=True, method=method ) - assert distance == 3 - assert qldpc.codes.get_distance_quantum(np.empty((0, 4), dtype=np.uint8), [[1, 0, 0, 0]]) == 3 - assert qldpc.codes.get_distance_quantum([[1, 0, 0, 0]], [[1, 0, 0, 0]]) == 3 + assert distance is no_distance + no_logical_ops = np.empty((0, 4), dtype=np.uint8) + assert qldpc.codes.get_distance_quantum(no_logical_ops, [[1, 0, 0, 0]]) is no_distance + assert qldpc.codes.get_distance_quantum([[1, 0, 0, 0]], [[1, 0, 0, 0]]) is no_distance - with pytest.raises(ValueError, match="incompatible widths"): + with pytest.raises(ValueError, match="The stabilizers have 2 columns"): qldpc.codes.get_distance_quantum([[1, 0, 0]], [[1, 0]], homogeneous=True) with pytest.raises(ValueError, match="even number of columns"): qldpc.codes.get_distance_quantum([[1, 0, 0]], [], homogeneous=False) @@ -187,10 +194,14 @@ def test_brouwer_zimmermann_special_inputs() -> None: np.empty((0, 3), dtype=np.uint8), homogeneous=False, ) - with pytest.raises(ValueError, match="two-dimensional"): + with pytest.raises(ValueError, match="The generators must be a two-dimensional"): qldpc.codes.get_distance_classical(np.zeros((1, 1, 1), dtype=np.uint8)) - with pytest.raises(ValueError, match="infer the width"): + with pytest.raises(ValueError, match="infer the width of generators"): qldpc.codes.get_distance_classical([]) + with pytest.raises(ValueError, match="The logical_ops must be a two-dimensional"): + qldpc.codes.get_distance_quantum(np.zeros((1, 1, 2), dtype=np.uint8), [[1, 0]]) + with pytest.raises(ValueError, match="The stabilizers must be a two-dimensional"): + qldpc.codes.get_distance_quantum([[1, 0]], np.zeros((1, 1, 2), dtype=np.uint8)) @pytest.mark.parametrize( @@ -343,7 +354,7 @@ def test_max_weight(request: pytest.FixtureRequest) -> None: ((word | word >> 4) & 0b1111).bit_count() for word in get_span(basis) - get_span(basis[:2]) ) - def check(distance_with_limit: Callable[..., int], distance: int, method: str) -> None: + def check(distance_with_limit: Callable[..., float], distance: int, method: str) -> None: assert distance >= 2 for max_weight in [distance - 2, distance]: result = distance_with_limit(cutoff=0, max_weight=max_weight, method=method) @@ -372,11 +383,18 @@ def test_max_weight_caps_cutoff() -> None: def test_max_weight_conventions_and_validation() -> None: - """The limit also caps the results returned for inputs without nontrivial words.""" - assert qldpc.codes.get_distance_classical([[0, 0, 0]], max_weight=1) == 2 + """The limit leaves NO_DISTANCE, for inputs without nontrivial words, unchanged.""" + no_distance = qldpc.codes.NO_DISTANCE + assert qldpc.codes.get_distance_classical([[0, 0, 0]], max_weight=1) is no_distance no_logical_ops = np.zeros((0, 4), dtype=int) - assert qldpc.codes.get_distance_quantum(no_logical_ops, [[1, 1, 0, 0]], max_weight=0) == 1 - assert qldpc.codes.get_distance_quantum([[1, 1, 0, 0]], [[1, 1, 0, 0]], max_weight=1) == 2 + assert ( + qldpc.codes.get_distance_quantum(no_logical_ops, [[1, 1, 0, 0]], max_weight=0) + is no_distance + ) + assert ( + qldpc.codes.get_distance_quantum([[1, 1, 0, 0]], [[1, 1, 0, 0]], max_weight=1) + is no_distance + ) with pytest.raises(ValueError, match="max_weight must be nonnegative"): qldpc.codes.get_distance_classical([[1, 1]], max_weight=-1) with pytest.raises(ValueError, match="max_weight must be nonnegative"): @@ -409,3 +427,20 @@ def test_distance_backend_validation() -> None: # The types of method names agree with the table that validates them. assert set(get_args(core.DistanceMethod)) == set(core._DISTANCE_METHODS) assert set(get_args(core._CanonicalDistanceMethod)) == set(core._DISTANCE_METHODS.values()) + + +def test_no_distance_behaves_as_nan_and_survives_copies() -> None: + """NO_DISTANCE is a float nan, and copies and pickles of it are the same member.""" + no_distance = core.NO_DISTANCE + assert qldpc.codes.NoDistance is qldpc.codes.distance.NoDistance is core.NoDistance + assert qldpc.codes.NO_DISTANCE is qldpc.codes.distance.NO_DISTANCE is core.NoDistance.NAN + assert isinstance(no_distance, float) and math.isnan(no_distance) + assert str(no_distance) == f"{no_distance}" == repr(no_distance) == "nan" + assert repr((2, 0, no_distance)) == "(2, 0, nan)" + assert not no_distance >= 3 and no_distance != core.NO_DISTANCE + restored = pickle.loads(pickle.dumps(no_distance)) # noqa: S301 - trusted in-memory round trip + assert restored is no_distance + assert copy.copy(no_distance) is copy.deepcopy(no_distance) is no_distance + assert core.NoDistance(math.nan) is no_distance + with pytest.raises(ValueError, match="not a valid NoDistance"): + core.NoDistance(1.0) diff --git a/src/qldpc/codes/distance/logical_counts.py b/src/qldpc/codes/distance/logical_counts.py index 4701b1c89..7b79123c8 100644 --- a/src/qldpc/codes/distance/logical_counts.py +++ b/src/qldpc/codes/distance/logical_counts.py @@ -45,6 +45,7 @@ from ._search.brouwer_zimmermann import get_distance_brouwer_zimmermann_many from ._search.brute_force import get_nested_distance_brute_force from ._tuning import DEFAULT_TUNING, SearchTuning, get_tuning +from .core import NO_DISTANCE, NoDistance # codes.common imports this package, and codes.quantum imports codes.common, so code classes are # imported only where they are used. @@ -78,26 +79,41 @@ # results and entry points -@dataclass(frozen=True) +@dataclass(frozen=True, eq=False) class CSSLogicalCounts: """Counts of the nontrivial logical operators of one weight in a binary CSS code. Attributes: weight: The weight of the counted operators. If no weight was requested, this is the - distance of the code, or the number of qubits if the code encodes no qubits, as with - :func:`qldpc.codes.distance.get_distance_quantum`. + distance of the code, or ``NO_DISTANCE`` (see :class:`~qldpc.codes.distance.NoDistance`) + if the code encodes no qubits, and so has no distance. count_x: The number of X-only nontrivial logical operators of weight ``weight``. count_z: The number of Z-only nontrivial logical operators of weight ``weight``. count_pauli: The number of nontrivial logical Pauli operators of weight ``weight``, ignoring phase, or None if it was not requested. It includes the X-only and Z-only operators counted by ``count_x`` and ``count_z``, as well as Y and mixed operators. + + Two counts of exactly the same class are equal if their attributes are, and ``NO_DISTANCE`` + equals itself here, although as a ``nan`` it does not equal itself elsewhere. """ - weight: int + weight: int | Literal[NoDistance.NAN] count_x: int count_z: int count_pauli: int | None + def __eq__(self, other: object) -> bool: + if type(other) is not type(self): + return NotImplemented + # Tuples compare identical items as equal, so NO_DISTANCE equals itself. + return self._as_tuple() == other._as_tuple() + + def __hash__(self) -> int: + return hash(self._as_tuple()) + + def _as_tuple(self) -> tuple[int | float | None, ...]: + return self.weight, self.count_x, self.count_z, self.count_pauli + class CSSLogicalCounter: """Exact counter of the logical operators of a given weight in a binary CSS code. @@ -223,6 +239,8 @@ def count_logical_ops( Args: weight: Weight of the operators to count, between zero and the number of qubits. If None (the default), count the operators whose weight is the distance of the code. + A code that encodes no qubits has no distance, and its counts are then zero at the + weight ``NO_DISTANCE``. method: ``"brouwer_zimmermann"`` (the default), ``"macwilliams"``, or ``"brute_force"``. MacWilliams uses check-space weight distributions, cached on the counter but not transferred through deformation. It does not support subsystem @@ -242,7 +260,7 @@ def count_logical_ops( Returns: The counted weight and the counts. If no weight is given and the code encodes no - qubits, the returned weight is the number of qubits and every count is zero. + qubits, the returned weight is ``NO_DISTANCE`` and every count is zero. Raises: TypeError: If the weight is not an integer, or is a bool. @@ -277,7 +295,7 @@ def count_logical_ops( if not self.dimension: empty_pauli = 0 if count_pauli else None return CSSLogicalCounts( - self.num_qubits if weight is None else int(weight), 0, 0, empty_pauli + NO_DISTANCE if weight is None else int(weight), 0, 0, empty_pauli ) if method == "macwilliams": return self._count_with_macwilliams( @@ -754,7 +772,7 @@ def _warm_up_compiled_functions() -> None: labels, basis_family, upper_bound, - *information_sets.get_overlapping_family_policy(DEFAULT_TUNING, len(basis)), + information_sets.get_overlapping_family_policy(DEFAULT_TUNING, len(basis)), ) diff --git a/src/qldpc/codes/distance/logical_counts_test.py b/src/qldpc/codes/distance/logical_counts_test.py index c21f33756..1c01e0c37 100644 --- a/src/qldpc/codes/distance/logical_counts_test.py +++ b/src/qldpc/codes/distance/logical_counts_test.py @@ -72,7 +72,7 @@ def test_known_codes(method: lc.LogicalCountMethod) -> None: # The Z sector is heavier, so it contributes nothing at the global minimum. ([[1, 1, 0]], [[0, 0, 1]], Counts(1, 2, 0, 2)), (np.zeros((0, 3)), np.zeros((0, 3)), Counts(1, 3, 3, 9)), - (np.eye(4), np.zeros((0, 4)), Counts(4, 0, 0, 0)), + (np.eye(4), np.zeros((0, 4)), Counts(codes.NO_DISTANCE, 0, 0, 0)), ) for checks_x, checks_z, expected in examples: counter = get_counter(checks_x, checks_z) @@ -205,6 +205,7 @@ def test_compile_search_setup(monkeypatch: pytest.MonkeyPatch) -> None: def get_results() -> list[qldpc.codes.CSSLogicalCounts]: counter = qldpc.codes.CSSLogicalCounter(code) results = [counter.count_logical_ops()] + assert results[0].weight is not qldpc.codes.NO_DISTANCE with use_tuning(max_information_sets=4): results.append(counter.count_logical_ops(count_pauli=True)) for deformed in (counter.deformed_cnot(0, 1), counter.deformed([0, 1], [[1, 0], [1, 1]])): @@ -335,6 +336,7 @@ def record(*args: Any, **kwargs: Any) -> information_sets.BasisFamily: # A count of the same weight, without the distance search, builds the families of the scans. expected = lc.CSSLogicalCounter(code).count_logical_ops() + assert expected.weight is not codes.NO_DISTANCE sizes_with_distance = collections.Counter(sizes) sizes.clear() lc.CSSLogicalCounter(code).count_logical_ops(expected.weight) @@ -616,7 +618,17 @@ def test_deformed_counters_walk_their_own_checks() -> None: def test_macwilliams_capabilities() -> None: """Unsupported codes and requests raise errors, and still work with the other methods.""" no_logicals = get_counter(np.eye(4)[:2], np.eye(4)[2:]) - assert no_logicals.count_logical_ops(method="macwilliams") == Counts(4, 0, 0, None) + no_distance = no_logicals.count_logical_ops(method="macwilliams") + assert no_distance == Counts(codes.NO_DISTANCE, 0, 0, None) + assert no_distance.weight is codes.NO_DISTANCE + assert hash(no_distance) == hash(Counts(codes.NO_DISTANCE, 0, 0, None)) + assert no_distance != Counts(4, 0, 0, None) and no_distance != (codes.NO_DISTANCE, 0, 0, None) + + # As for dataclasses, counts of different classes differ. + class SubclassCounts(Counts): + pass + + assert no_distance != SubclassCounts(codes.NO_DISTANCE, 0, 0, None) extended_hamming = codes.ExtendedHammingCode(3).matrix self_dual = get_counter(extended_hamming, extended_hamming) assert self_dual.dimension == 0 diff --git a/src/qldpc/codes/quantum/hypergraph_product.py b/src/qldpc/codes/quantum/hypergraph_product.py index d4e00b2dc..58cc9f2ea 100644 --- a/src/qldpc/codes/quantum/hypergraph_product.py +++ b/src/qldpc/codes/quantum/hypergraph_product.py @@ -7,7 +7,7 @@ import collections import itertools from collections.abc import Collection, Iterable, Sequence -from typing import TypeVar +from typing import Literal, TypeVar import galois import numpy as np @@ -35,7 +35,7 @@ SimplexCode, ) from ..common import ClassicalCode, CSSCode -from ..distance import DistanceMethod +from ..distance import NO_DISTANCE, DistanceMethod, NoDistance #################################################################################################### # hypergraph product codes and their variants @@ -378,7 +378,9 @@ def get_canonical_logical_line_ops( logical_ops_z = scipy.linalg.block_diag(logical_ops_z_l, logical_ops_z_r) return logical_ops_x.view(field), logical_ops_z.view(field) - def _get_distance_exact(self, pauli: PauliXZ | None, *, method: DistanceMethod) -> int | float: + def _get_distance_exact( + self, pauli: PauliXZ | None, *, method: DistanceMethod + ) -> int | Literal[NoDistance.NAN]: """Exact distance calculation for hypergraph product codes. These calculations are based on arXiv:2308.15520, but additionally allow for the separate @@ -401,15 +403,20 @@ def _get_distance_exact(self, pauli: PauliXZ | None, *, method: DistanceMethod) if pauli is Pauli.X: transpose_a = ClassicalCode(self.code_a.matrix.T) - dist_a = transpose_a.get_distance(method=method) if dependent_b else np.nan - dist_b = self.code_b.get_distance(method=method) if self.code_a.dimension else np.nan + dist_a = transpose_a.get_distance(method=method) if dependent_b else NO_DISTANCE + dist_b = ( + self.code_b.get_distance(method=method) if self.code_a.dimension else NO_DISTANCE + ) else: assert pauli is Pauli.Z - dist_a = self.code_a.get_distance(method=method) if self.code_b.dimension else np.nan + dist_a = ( + self.code_a.get_distance(method=method) if self.code_b.dimension else NO_DISTANCE + ) transpose_b = ClassicalCode(self.code_b.matrix.T) - dist_b = transpose_b.get_distance(method=method) if dependent_a else np.nan + dist_b = transpose_b.get_distance(method=method) if dependent_a else NO_DISTANCE - return dist_a if np.isnan(dist_b) else dist_b if np.isnan(dist_a) else min(dist_a, dist_b) + distances = [distance for distance in (dist_a, dist_b) if distance is not NO_DISTANCE] + return min(distances, default=NO_DISTANCE) class CHGPCode(HGPCode): @@ -462,7 +469,7 @@ def __init__( """ cyclic_code = CyclicCode(bits, poly, field) distance = cyclic_code.get_distance() - if np.isnan(distance): + if distance is NO_DISTANCE: raise ValueError( "Cannot construct a CRCode from a cyclic code without a defined distance" ) @@ -558,7 +565,9 @@ def get_canonical_logical_line_ops( logical_ops_z = np.kron(generator_a, pivots_b) return logical_ops_x.view(field), logical_ops_z.view(field) - def _get_distance_exact(self, pauli: PauliXZ | None, *, method: DistanceMethod) -> int | float: + def _get_distance_exact( + self, pauli: PauliXZ | None, *, method: DistanceMethod + ) -> int | Literal[NoDistance.NAN]: """Exact distance calculation for subsystem hypergraph product codes.""" match pauli: case Pauli.X: diff --git a/src/qldpc/codes/quantum/self_orthogonal.py b/src/qldpc/codes/quantum/self_orthogonal.py index bb07987bf..18ca10c79 100644 --- a/src/qldpc/codes/quantum/self_orthogonal.py +++ b/src/qldpc/codes/quantum/self_orthogonal.py @@ -4,6 +4,8 @@ from __future__ import annotations +from typing import Literal + import galois import numpy as np import scipy @@ -17,7 +19,7 @@ ReedMullerCode, ) from ..common import CSSCode -from ..distance import DistanceMethod +from ..distance import DistanceMethod, NoDistance #################################################################################################### # self-orthogonal CSS codes @@ -127,7 +129,9 @@ def _assert_valid_params(order: int, size: int) -> None: + f"Provided: (r, m) = ({order}, {size})" ) - def _get_distance_exact(self, pauli: PauliXZ | None, *, method: DistanceMethod) -> int | float: + def _get_distance_exact( + self, pauli: PauliXZ | None, *, method: DistanceMethod + ) -> int | Literal[NoDistance.NAN]: """Closed-form distance of CSS(RM(r, m), RM(r, m)) codes: 2**(r + 1).""" return 2 ** (self._order + 1) diff --git a/src/qldpc/external/_codes.py b/src/qldpc/external/_codes.py index 790facf48..cde7e2b47 100644 --- a/src/qldpc/external/_codes.py +++ b/src/qldpc/external/_codes.py @@ -239,7 +239,9 @@ def get_qdistrnd_distance_bound( bound = "".join([line for line in lines if not line.startswith("#")]) if not bound: raise ValueError(f"Could not parse a distance bound from QDistRnd output: {output!r}") - return int(bound) + + # QDistRnd negates a bound at or below its argument mindist, the cutoff, to mark an early stop. + return abs(int(bound)) def _get_gap_sparse_matrix_commands( diff --git a/src/qldpc/external/_codes_test.py b/src/qldpc/external/_codes_test.py index b264d47ec..2b70a6c57 100644 --- a/src/qldpc/external/_codes_test.py +++ b/src/qldpc/external/_codes_test.py @@ -133,6 +133,10 @@ def test_get_qdistrnd_distance_bound() -> None: == 3 ) + # QDistRnd negates a bound at or below the cutoff when it stops early + with unittest.mock.patch("qldpc.external.gap.get_output", return_value="-2"): + assert _codes.get_qdistrnd_distance_bound(matrix, cutoff=2) == 2 + # QDistRnd produced no output at all with ( unittest.mock.patch("qldpc.external.gap.get_output", return_value=""),