From f2e1112e99c1ad3d08b05a9ae2b00c2e9445c4a8 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Fri, 6 Dec 2024 11:28:08 -0800 Subject: [PATCH 01/21] Take a first pass at bitpacked fields. --- src/galois/_domains/_array.py | 22 ++++- src/galois/_domains/_function.py | 7 +- src/galois/_domains/_meta.py | 1 + src/galois/_fields/_gf2.py | 143 ++++++++++++++++++++++++++++++- 4 files changed, 169 insertions(+), 4 deletions(-) diff --git a/src/galois/_domains/_array.py b/src/galois/_domains/_array.py index 283de9d49..1a5155d01 100644 --- a/src/galois/_domains/_array.py +++ b/src/galois/_domains/_array.py @@ -53,10 +53,15 @@ def __new__( dtype = cls._get_dtype(dtype) x = cls._verify_array_like_types_and_values(x) - array = np.array(x, dtype=dtype, copy=copy, order=order, ndmin=ndmin) + if cls._bitpacked: + array = np.packbits(x, axis=-1) + else: + array = np.array(x, dtype=dtype, copy=copy, order=order, ndmin=ndmin) # Perform view without verification since the elements were verified in _verify_array_like_types_and_values() - return cls._view(array) + galois_array = cls._view(array) + galois_array.original_shape = x.shape + return galois_array @classmethod def _get_dtype(cls, dtype: DTypeLike) -> np.dtype: @@ -172,6 +177,9 @@ def Zeros(cls, shape: ShapeLike, dtype: DTypeLike | None = None) -> Self: """ dtype = cls._get_dtype(dtype) array = np.zeros(shape, dtype=dtype) + if cls._bitpacked: + array = np.packbits(array, axis=-1) + return cls._view(array) @classmethod @@ -190,6 +198,9 @@ def Ones(cls, shape: ShapeLike, dtype: DTypeLike | None = None) -> Self: """ dtype = cls._get_dtype(dtype) array = np.ones(shape, dtype=dtype) + if cls._bitpacked: + array = np.packbits(array, axis=-1) + return cls._view(array) @classmethod @@ -227,6 +238,8 @@ def Range( raise ValueError(f"Argument 'stop' must be within the field's order {cls.order}, not {stop}.") array = np.arange(start, stop, step=step, dtype=dtype) + if cls._bitpacked: + array = np.packbits(array, axis=-1) return cls._view(array) @@ -295,6 +308,8 @@ def Random( for _ in iterator: array[iterator.multi_index] = random.randint(low, high - 1) + if cls._bitpacked: + array = np.packbits(array, axis=-1) return cls._view(array) @classmethod @@ -313,6 +328,9 @@ def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: """ dtype = cls._get_dtype(dtype) array = np.identity(size, dtype=dtype) + if cls._bitpacked: + array = np.packbits(array, axis=-1) + return cls._view(array) ############################################################################### diff --git a/src/galois/_domains/_function.py b/src/galois/_domains/_function.py index eaf7d4a9f..d1cd9c807 100644 --- a/src/galois/_domains/_function.py +++ b/src/galois/_domains/_function.py @@ -334,6 +334,7 @@ class FunctionMixin(np.ndarray, metaclass=ArrayMeta): np.convolve: "_convolve", np.fft.fft: "_fft", np.fft.ifft: "_ifft", + np.array_equal: "_array_equal", } _convolve: Function @@ -353,7 +354,10 @@ def __array_function__(self, func, types, args, kwargs): field = type(self) if func in field._OVERRIDDEN_FUNCTIONS: - output = getattr(field, field._OVERRIDDEN_FUNCTIONS[func])(*args, **kwargs) + try: + output = getattr(field, field._OVERRIDDEN_FUNCTIONS[func])(*args, **kwargs) + except AttributeError: + output = super().__array_function__(func, types, args, kwargs) elif func in field._UNSUPPORTED_FUNCTIONS: raise NotImplementedError( @@ -377,6 +381,7 @@ def __array_function__(self, func, types, args, kwargs): return output + def dot(self, b, out=None): # The `np.dot(a, b)` ufunc is also available as `a.dot(b)`. Need to override this method for # consistent results. diff --git a/src/galois/_domains/_meta.py b/src/galois/_domains/_meta.py index 892bf549f..a15557416 100644 --- a/src/galois/_domains/_meta.py +++ b/src/galois/_domains/_meta.py @@ -35,6 +35,7 @@ def __init__(cls, name, bases, namespace, **kwargs): cls._irreducible_poly_int: int = kwargs.get("irreducible_poly_int", 0) cls._primitive_element: int = kwargs.get("primitive_element", 0) cls._dtypes = cls._determine_dtypes() + cls._bitpacked = kwargs.get("bitpacked", False) if cls._dtypes == [np.object_]: cls._default_ufunc_mode = "python-calculate" diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 212feed11..ccc916fe3 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -17,7 +17,7 @@ sqrt_ufunc, subtract_ufunc, ) -from .._domains._ufunc import UFuncMixin +from .._domains._ufunc import UFuncMixin, matmul_ufunc from .._helper import export from ._array import FieldArray @@ -100,6 +100,98 @@ def __init_subclass__(cls) -> None: cls._log = log(cls) cls._sqrt = sqrt(cls) +class add_ufunc_bitpacked(add_ufunc): + """ + Addition ufunc dispatcher w/ support for bit-packed fields. + """ + def __call__(self, ufunc, method, inputs, kwargs, meta): + output = super().__call__(ufunc, method, inputs, kwargs, meta) + output.original_shape = inputs[0].original_shape + return output + +class subtract_ufunc_bitpacked(subtract_ufunc): + """ + Subtraction ufunc dispatcher w/ support for bit-packed fields. + """ + def __call__(self, ufunc, method, inputs, kwargs, meta): + output = super().__call__(ufunc, method, inputs, kwargs, meta) + output.original_shape = inputs[0].original_shape + return output + +class multiply_ufunc_bitpacked(multiply_ufunc): + """ + Multiply ufunc dispatcher w/ support for bit-packed fields. + """ + def __call__(self, ufunc, method, inputs, kwargs, meta): + output = super().__call__(ufunc, method, inputs, kwargs, meta) + output.original_shape = inputs[0].original_shape + return output + +class divide_ufunc_bitpacked(divide): + """ + Divide ufunc dispatcher w/ support for bit-packed fields. + """ + def __call__(self, ufunc, method, inputs, kwargs, meta): + output = super().__call__(ufunc, method, inputs, kwargs, meta) + output.original_shape = inputs[0].original_shape + return output + +class matmul_ufunc_bitpacked(matmul_ufunc): + """ + Matmul ufunc dispatcher w/ support for bit-packed fields. + """ + def __call__(self, ufunc, method, inputs, kwargs, meta): + a, b = inputs + + if hasattr(a, "original_shape"): + a = np.unpackbits(a.view(np.ndarray), axis=-1, count=a.original_shape[-1]).view(GF2BP) + if hasattr(b, "original_shape"): + b = np.unpackbits(b.view(np.ndarray), axis=-1, count=b.original_shape[-1]).view(GF2BP) + + inputs = (a, b) + output = super().__call__(ufunc, method, inputs, kwargs, meta) + original_shape = output.shape + output = self.field._view(np.packbits(output.view(np.ndarray), axis=-1)) + output.original_shape = original_shape + + return output + + +def array_equal_bitpacked(a: FieldArray, b: FieldArray) -> bool: + if a.shape != b.shape: + if hasattr(a, "original_shape"): + a = np.unpackbits(a.view(np.ndarray), axis=-1, count=a.original_shape[-1]) + + if hasattr(b, "original_shape"): + b = np.unpackbits(b.view(np.ndarray), axis=-1, count=b.original_shape[-1]) + + return np.core.numeric.array_equal(a, b) + + +class UFuncMixin_2_1_BitPacked(UFuncMixin): + """ + A mixin class that provides explicit calculation arithmetic for GF(2). + """ + + def __init_subclass__(cls) -> None: + super().__init_subclass__() + cls._add = add_ufunc_bitpacked(cls, override=np.bitwise_xor) + cls._negative = negative_ufunc(cls, override=np.positive) + cls._subtract = subtract_ufunc_bitpacked(cls, override=np.bitwise_xor) + cls._multiply = multiply_ufunc_bitpacked(cls, override=np.bitwise_and) + cls._reciprocal = reciprocal(cls) + cls._divide = divide_ufunc_bitpacked(cls) + cls._power = power(cls) + cls._log = log(cls) + cls._sqrt = sqrt(cls) + cls._array_equal = array_equal_bitpacked + + @classmethod + def _assign_ufuncs(cls): + super()._assign_ufuncs() + + # We have to set this here because ArrayMeta would override it. + cls._matmul = matmul_ufunc_bitpacked(cls) # NOTE: There is a "verbatim" block in the docstring because we were not able to monkey-patch GF2 like the # other classes in docs/conf.py. So, technically, at doc-build-time issubclass(galois.GF2, galois.FieldArray) == False @@ -116,7 +208,52 @@ class GF2( order=2, irreducible_poly_int=3, is_primitive_poly=True, + primitive_element=1 +): + r""" + A :obj:`~galois.FieldArray` subclass over $\mathrm{GF}(2)$. + + .. info:: + + This class is a pre-generated :obj:`~galois.FieldArray` subclass generated with `galois.GF(2)` and is + included in the API for convenience. + + Examples: + This class is equivalent, and in fact identical, to the :obj:`~galois.FieldArray` subclass returned from the + class factory :func:`~galois.GF`. + + .. ipython:: + + In [2]: galois.GF2 is galois.GF(2) + + @verbatim + In [3]: issubclass(galois.GF2, galois.FieldArray) + Out[3]: True + + In [4]: print(galois.GF2.properties) + + Create a :obj:`~galois.FieldArray` instance using :obj:`~galois.GF2`'s constructor. + + .. ipython:: python + + x = galois.GF2([1, 0, 1, 1]); x + isinstance(x, galois.GF2) + + Group: + galois-fields + """ + +@export +class GF2BP( + FieldArray, + UFuncMixin_2_1_BitPacked, + characteristic=2, + degree=1, + order=2, + irreducible_poly_int=3, + is_primitive_poly=True, primitive_element=1, + bitpacked=True, ): r""" A :obj:`~galois.FieldArray` subclass over $\mathrm{GF}(2)$. @@ -155,3 +292,7 @@ class factory :func:`~galois.GF`. GF2._default_ufunc_mode = "jit-calculate" GF2._ufunc_modes = ["jit-calculate", "python-calculate"] GF2.compile("auto") + +GF2BP._default_ufunc_mode = "jit-calculate" +GF2BP._ufunc_modes = ["jit-calculate", "python-calculate"] +GF2BP.compile("auto") From a209f8630a4d7100a2f1b49443ae331b8c1ccafb Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Tue, 10 Dec 2024 13:38:49 -0800 Subject: [PATCH 02/21] Check in WIP on implicit bitpacking --- src/galois/_domains/_array.py | 51 +++++++++++++++++++++++++++----- src/galois/_domains/_function.py | 1 + src/galois/_domains/_meta.py | 2 ++ src/galois/_fields/_array.py | 40 +++++++++++++++++++++++++ src/galois/_fields/_gf2.py | 38 ++++++++++++++++++++---- 5 files changed, 120 insertions(+), 12 deletions(-) diff --git a/src/galois/_domains/_array.py b/src/galois/_domains/_array.py index 1a5155d01..8aaef916b 100644 --- a/src/galois/_domains/_array.py +++ b/src/galois/_domains/_array.py @@ -53,14 +53,18 @@ def __new__( dtype = cls._get_dtype(dtype) x = cls._verify_array_like_types_and_values(x) - if cls._bitpacked: + packed = False + if cls._bitpacked and x.size > 0: array = np.packbits(x, axis=-1) + packed = True else: array = np.array(x, dtype=dtype, copy=copy, order=order, ndmin=ndmin) # Perform view without verification since the elements were verified in _verify_array_like_types_and_values() galois_array = cls._view(array) - galois_array.original_shape = x.shape + if cls._bitpacked: + galois_array.original_shape = x.shape + galois_array.packed = packed return galois_array @classmethod @@ -178,9 +182,16 @@ def Zeros(cls, shape: ShapeLike, dtype: DTypeLike | None = None) -> Self: dtype = cls._get_dtype(dtype) array = np.zeros(shape, dtype=dtype) if cls._bitpacked: + original_shape = array.shape array = np.packbits(array, axis=-1) - return cls._view(array) + array = cls._view(array) + + if cls._bitpacked: + array.packed = True + array.original_shape = original_shape + + return array @classmethod def Ones(cls, shape: ShapeLike, dtype: DTypeLike | None = None) -> Self: @@ -199,10 +210,16 @@ def Ones(cls, shape: ShapeLike, dtype: DTypeLike | None = None) -> Self: dtype = cls._get_dtype(dtype) array = np.ones(shape, dtype=dtype) if cls._bitpacked: + original_shape = array.shape array = np.packbits(array, axis=-1) - return cls._view(array) + array = cls._view(array) + + if cls._bitpacked: + array.packed = True + array.original_shape = original_shape + return array @classmethod def Range( cls, @@ -239,10 +256,16 @@ def Range( array = np.arange(start, stop, step=step, dtype=dtype) if cls._bitpacked: + original_shape = array.shape array = np.packbits(array, axis=-1) - return cls._view(array) + array = cls._view(array) + + if cls._bitpacked: + array.packed = True + array.original_shape = original_shape + return array @classmethod def Random( cls, @@ -309,9 +332,16 @@ def Random( array[iterator.multi_index] = random.randint(low, high - 1) if cls._bitpacked: + original_shape = array.shape array = np.packbits(array, axis=-1) - return cls._view(array) + array = cls._view(array) + + if cls._bitpacked: + array.packed = True + array.original_shape = original_shape + + return array @classmethod def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: r""" @@ -329,9 +359,16 @@ def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: dtype = cls._get_dtype(dtype) array = np.identity(size, dtype=dtype) if cls._bitpacked: + original_shape = array.shape array = np.packbits(array, axis=-1) - return cls._view(array) + array = cls._view(array) + + if cls._bitpacked: + array.packed = True + array.original_shape = original_shape + + return array ############################################################################### # Ufunc compilation routines diff --git a/src/galois/_domains/_function.py b/src/galois/_domains/_function.py index d1cd9c807..c7fdf2cd6 100644 --- a/src/galois/_domains/_function.py +++ b/src/galois/_domains/_function.py @@ -335,6 +335,7 @@ class FunctionMixin(np.ndarray, metaclass=ArrayMeta): np.fft.fft: "_fft", np.fft.ifft: "_ifft", np.array_equal: "_array_equal", + np.concatenate: "_concatenate", } _convolve: Function diff --git a/src/galois/_domains/_meta.py b/src/galois/_domains/_meta.py index a15557416..85d726975 100644 --- a/src/galois/_domains/_meta.py +++ b/src/galois/_domains/_meta.py @@ -36,6 +36,8 @@ def __init__(cls, name, bases, namespace, **kwargs): cls._primitive_element: int = kwargs.get("primitive_element", 0) cls._dtypes = cls._determine_dtypes() cls._bitpacked = kwargs.get("bitpacked", False) + if cls._bitpacked and (cls._order != 2 or cls._degree != 1): + raise ValueError("Only fields over GF(2) support bit-packing.") if cls._dtypes == [np.object_]: cls._default_ufunc_mode = "python-calculate" diff --git a/src/galois/_fields/_array.py b/src/galois/_fields/_array.py index 3a074bc92..643aad463 100644 --- a/src/galois/_fields/_array.py +++ b/src/galois/_fields/_array.py @@ -7,6 +7,7 @@ from typing import Generator import numpy as np +import numpy.typing as npt from typing_extensions import Literal, Self from .._domains import Array, _linalg @@ -960,6 +961,45 @@ def primitive_roots_of_unity(cls, n: int) -> Self: # Instance methods ############################################################################### + @property + def shape(self) -> npt._shape: + if self._bitpacked and hasattr(self, "original_shape"): + return self.original_shape + + return super().shape + + @shape.setter + def shape(self, value: npt._shape) -> None: + if self._bitpacked: + raise NotImplementedError("You cannot set the shape on a bit-packed array.") + + super().shape = value + + def __array__(self) -> np.ndarray: + if self._bitpacked and hasattr(self, "original_shape"): + return np.unpackbits(self.view(np.ndarray), axis=-1, count=self.original_shape[-1]) + + return self + + def __getitem__(self, index): + return self.__array__()[index] + + def __setitem__(self, index, value): + array_unpacked = self.__array__() + array_unpacked[index] = value + self.view(np.ndarray)[:] = np.packbits(array_unpacked, axis=-1)[:] + + @property + def T(self) -> Self: + if self._bitpacked and hasattr(self, "original_shape"): + transposed = np.unpackbits(self.view(np.ndarray), axis=-1, count=self.original_shape[-1]).T + original_shape = transposed.shape + transposed = np.packbits(transposed, axis=-1) + transposed.original_shape = original_shape + return transposed + + return super().T + def additive_order(self) -> int | np.ndarray: r""" Computes the additive order of each element in $x$. diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index ccc916fe3..8f30f8185 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -145,8 +145,12 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): if hasattr(a, "original_shape"): a = np.unpackbits(a.view(np.ndarray), axis=-1, count=a.original_shape[-1]).view(GF2BP) + else: + a = a.view(GF2BP) if hasattr(b, "original_shape"): b = np.unpackbits(b.view(np.ndarray), axis=-1, count=b.original_shape[-1]).view(GF2BP) + else: + b = b.view(GF2BP) inputs = (a, b) output = super().__call__(ufunc, method, inputs, kwargs, meta) @@ -158,15 +162,37 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): def array_equal_bitpacked(a: FieldArray, b: FieldArray) -> bool: - if a.shape != b.shape: - if hasattr(a, "original_shape"): - a = np.unpackbits(a.view(np.ndarray), axis=-1, count=a.original_shape[-1]) + unpack_a = False + unpack_b = False - if hasattr(b, "original_shape"): - b = np.unpackbits(b.view(np.ndarray), axis=-1, count=b.original_shape[-1]) + a_is_bitpacked = hasattr(a, "original_shape") + b_is_bitpacked = hasattr(b, "original_shape") + if a_is_bitpacked and b_is_bitpacked and a.shape != b.shape: + unpack_a = True + unpack_b = True + elif a_is_bitpacked: + unpack_a = True + elif b_is_bitpacked: + unpack_b = True + + if unpack_a: + a = np.unpackbits(a.view(np.ndarray), axis=-1, count=a.original_shape[-1]).view(GF2) + + if unpack_b: + b = np.unpackbits(b.view(np.ndarray), axis=-1, count=b.original_shape[-1]).view(GF2) return np.core.numeric.array_equal(a, b) +def concatenate_bitpacked(arrays, axis=None, out=None, **kwargs): + array_list = list(arrays) + for i, array in enumerate(arrays): + if hasattr(array, "original_shape"): + array_list[i] = np.unpackbits(array.view(np.ndarray), axis=-1, count=array.original_shape[-1]).view(np.ndarray) + else: + array_list[i] = array.view(np.ndarray) + + return np.core.multiarray.concatenate(tuple(array_list), axis=axis, out=out, **kwargs) + class UFuncMixin_2_1_BitPacked(UFuncMixin): """ @@ -185,6 +211,8 @@ def __init_subclass__(cls) -> None: cls._log = log(cls) cls._sqrt = sqrt(cls) cls._array_equal = array_equal_bitpacked + cls._concatenate = concatenate_bitpacked + cls._inv = None @classmethod def _assign_ufuncs(cls): From 3ec5def904cb63cd67d7888ec3f71b9ee85a3ea5 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Tue, 10 Dec 2024 16:33:38 -0800 Subject: [PATCH 03/21] Clean up formatting and move towards a more explicit version. --- src/galois/_domains/_array.py | 24 ++++--- src/galois/_domains/_function.py | 2 - src/galois/_fields/_array.py | 36 +++++----- src/galois/_fields/_gf2.py | 113 +++++++++++++++++-------------- 4 files changed, 95 insertions(+), 80 deletions(-) diff --git a/src/galois/_domains/_array.py b/src/galois/_domains/_array.py index 8aaef916b..721eacb35 100644 --- a/src/galois/_domains/_array.py +++ b/src/galois/_domains/_array.py @@ -7,6 +7,7 @@ import abc import contextlib import random +from numbers import Integral from typing import Generator import numpy as np @@ -54,17 +55,17 @@ def __new__( x = cls._verify_array_like_types_and_values(x) packed = False - if cls._bitpacked and x.size > 0: - array = np.packbits(x, axis=-1) + if cls._bitpacked and isinstance(x, np.ndarray) and x.size > 0: + array = np.packbits(x.view(np.ndarray), axis=-1) packed = True else: array = np.array(x, dtype=dtype, copy=copy, order=order, ndmin=ndmin) # Perform view without verification since the elements were verified in _verify_array_like_types_and_values() galois_array = cls._view(array) - if cls._bitpacked: - galois_array.original_shape = x.shape - galois_array.packed = packed + if cls._bitpacked and packed: + galois_array._unpacked_shape = x.shape + galois_array.packed = packed return galois_array @classmethod @@ -189,7 +190,7 @@ def Zeros(cls, shape: ShapeLike, dtype: DTypeLike | None = None) -> Self: if cls._bitpacked: array.packed = True - array.original_shape = original_shape + array._unpacked_shape = original_shape return array @@ -217,9 +218,10 @@ def Ones(cls, shape: ShapeLike, dtype: DTypeLike | None = None) -> Self: if cls._bitpacked: array.packed = True - array.original_shape = original_shape + array._unpacked_shape = original_shape return array + @classmethod def Range( cls, @@ -263,9 +265,10 @@ def Range( if cls._bitpacked: array.packed = True - array.original_shape = original_shape + array._unpacked_shape = original_shape return array + @classmethod def Random( cls, @@ -339,9 +342,10 @@ def Random( if cls._bitpacked: array.packed = True - array.original_shape = original_shape + array._unpacked_shape = original_shape return array + @classmethod def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: r""" @@ -366,7 +370,7 @@ def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: if cls._bitpacked: array.packed = True - array.original_shape = original_shape + array._unpacked_shape = original_shape return array diff --git a/src/galois/_domains/_function.py b/src/galois/_domains/_function.py index c7fdf2cd6..30c98d576 100644 --- a/src/galois/_domains/_function.py +++ b/src/galois/_domains/_function.py @@ -335,7 +335,6 @@ class FunctionMixin(np.ndarray, metaclass=ArrayMeta): np.fft.fft: "_fft", np.fft.ifft: "_ifft", np.array_equal: "_array_equal", - np.concatenate: "_concatenate", } _convolve: Function @@ -382,7 +381,6 @@ def __array_function__(self, func, types, args, kwargs): return output - def dot(self, b, out=None): # The `np.dot(a, b)` ufunc is also available as `a.dot(b)`. Need to override this method for # consistent results. diff --git a/src/galois/_fields/_array.py b/src/galois/_fields/_array.py index 643aad463..c995039dd 100644 --- a/src/galois/_fields/_array.py +++ b/src/galois/_fields/_array.py @@ -963,8 +963,8 @@ def primitive_roots_of_unity(cls, n: int) -> Self: @property def shape(self) -> npt._shape: - if self._bitpacked and hasattr(self, "original_shape"): - return self.original_shape + if self._bitpacked and hasattr(self, "_unpacked_shape"): + return self._unpacked_shape return super().shape @@ -975,24 +975,26 @@ def shape(self, value: npt._shape) -> None: super().shape = value - def __array__(self) -> np.ndarray: - if self._bitpacked and hasattr(self, "original_shape"): - return np.unpackbits(self.view(np.ndarray), axis=-1, count=self.original_shape[-1]) - - return self - - def __getitem__(self, index): - return self.__array__()[index] - - def __setitem__(self, index, value): - array_unpacked = self.__array__() - array_unpacked[index] = value - self.view(np.ndarray)[:] = np.packbits(array_unpacked, axis=-1)[:] + # def __array__(self) -> np.ndarray: + # if self._bitpacked and hasattr(self, "_unpacked_shape"): + # return np.unpackbits(self.view(np.ndarray), axis=-1, count=self._unpacked_shape[-1]) + # + # return self + # + # def __getitem__(self, index): + # return self.__array__()[index] + # + # def __setitem__(self, index, value): + # array_unpacked = self.__array__() + # array_unpacked[index] = value + # self.view(np.ndarray)[:] = np.packbits(array_unpacked, axis=-1)[:] @property def T(self) -> Self: - if self._bitpacked and hasattr(self, "original_shape"): - transposed = np.unpackbits(self.view(np.ndarray), axis=-1, count=self.original_shape[-1]).T + if self._bitpacked and hasattr(self, "_unpacked_shape"): + transposed = np.unpackbits( + self.view(np.ndarray), axis=-1, count=self._unpacked_shape[-1] + ).T original_shape = transposed.shape transposed = np.packbits(transposed, axis=-1) transposed.original_shape = original_shape diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 8f30f8185..78c49a036 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -4,8 +4,10 @@ from __future__ import annotations +import numpy import numpy as np +from ..typing import ElementLike, ArrayLike, DTypeLike from .._domains._lookup import ( add_ufunc, divide_ufunc, @@ -100,98 +102,79 @@ def __init_subclass__(cls) -> None: cls._log = log(cls) cls._sqrt = sqrt(cls) + class add_ufunc_bitpacked(add_ufunc): """ Addition ufunc dispatcher w/ support for bit-packed fields. """ + def __call__(self, ufunc, method, inputs, kwargs, meta): output = super().__call__(ufunc, method, inputs, kwargs, meta) - output.original_shape = inputs[0].original_shape + output._unpacked_shape = inputs[0]._unpacked_shape return output + class subtract_ufunc_bitpacked(subtract_ufunc): """ Subtraction ufunc dispatcher w/ support for bit-packed fields. """ + def __call__(self, ufunc, method, inputs, kwargs, meta): output = super().__call__(ufunc, method, inputs, kwargs, meta) - output.original_shape = inputs[0].original_shape + output._unpacked_shape = inputs[0]._unpacked_shape return output + class multiply_ufunc_bitpacked(multiply_ufunc): """ Multiply ufunc dispatcher w/ support for bit-packed fields. """ + def __call__(self, ufunc, method, inputs, kwargs, meta): output = super().__call__(ufunc, method, inputs, kwargs, meta) - output.original_shape = inputs[0].original_shape + output._unpacked_shape = inputs[0]._unpacked_shape return output + class divide_ufunc_bitpacked(divide): """ Divide ufunc dispatcher w/ support for bit-packed fields. """ + def __call__(self, ufunc, method, inputs, kwargs, meta): output = super().__call__(ufunc, method, inputs, kwargs, meta) - output.original_shape = inputs[0].original_shape + output._unpacked_shape = inputs[0]._unpacked_shape return output + class matmul_ufunc_bitpacked(matmul_ufunc): """ Matmul ufunc dispatcher w/ support for bit-packed fields. """ + def __call__(self, ufunc, method, inputs, kwargs, meta): a, b = inputs - if hasattr(a, "original_shape"): - a = np.unpackbits(a.view(np.ndarray), axis=-1, count=a.original_shape[-1]).view(GF2BP) - else: - a = a.view(GF2BP) - if hasattr(b, "original_shape"): - b = np.unpackbits(b.view(np.ndarray), axis=-1, count=b.original_shape[-1]).view(GF2BP) - else: - b = b.view(GF2BP) + assert isinstance(a, GF2BP) and isinstance(b, GF2BP) + assert a.shape[-1] == b.shape[0] + + # a has rows packed, so repack b's columns + b = np.packbits( + np.unpackbits(b.view(np.ndarray), axis=-1, count=b._unpacked_shape[-1]), + axis=0, + ).view(GF2BP) inputs = (a, b) output = super().__call__(ufunc, method, inputs, kwargs, meta) - original_shape = output.shape + unpacked_shape = output.shape output = self.field._view(np.packbits(output.view(np.ndarray), axis=-1)) - output.original_shape = original_shape + output._unpacked_shape = unpacked_shape return output -def array_equal_bitpacked(a: FieldArray, b: FieldArray) -> bool: - unpack_a = False - unpack_b = False - - a_is_bitpacked = hasattr(a, "original_shape") - b_is_bitpacked = hasattr(b, "original_shape") - if a_is_bitpacked and b_is_bitpacked and a.shape != b.shape: - unpack_a = True - unpack_b = True - elif a_is_bitpacked: - unpack_a = True - elif b_is_bitpacked: - unpack_b = True - - if unpack_a: - a = np.unpackbits(a.view(np.ndarray), axis=-1, count=a.original_shape[-1]).view(GF2) - - if unpack_b: - b = np.unpackbits(b.view(np.ndarray), axis=-1, count=b.original_shape[-1]).view(GF2) - - return np.core.numeric.array_equal(a, b) - -def concatenate_bitpacked(arrays, axis=None, out=None, **kwargs): - array_list = list(arrays) - for i, array in enumerate(arrays): - if hasattr(array, "original_shape"): - array_list[i] = np.unpackbits(array.view(np.ndarray), axis=-1, count=array.original_shape[-1]).view(np.ndarray) - else: - array_list[i] = array.view(np.ndarray) - - return np.core.multiarray.concatenate(tuple(array_list), axis=axis, out=out, **kwargs) +def not_implemented(*args, **kwargs): + return NotImplemented class UFuncMixin_2_1_BitPacked(UFuncMixin): @@ -210,16 +193,14 @@ def __init_subclass__(cls) -> None: cls._power = power(cls) cls._log = log(cls) cls._sqrt = sqrt(cls) - cls._array_equal = array_equal_bitpacked - cls._concatenate = concatenate_bitpacked - cls._inv = None @classmethod def _assign_ufuncs(cls): super()._assign_ufuncs() # We have to set this here because ArrayMeta would override it. - cls._matmul = matmul_ufunc_bitpacked(cls) + cls._matmul = matmul_ufunc_bitpacked(cls) + # NOTE: There is a "verbatim" block in the docstring because we were not able to monkey-patch GF2 like the # other classes in docs/conf.py. So, technically, at doc-build-time issubclass(galois.GF2, galois.FieldArray) == False @@ -236,7 +217,7 @@ class GF2( order=2, irreducible_poly_int=3, is_primitive_poly=True, - primitive_element=1 + primitive_element=1, ): r""" A :obj:`~galois.FieldArray` subclass over $\mathrm{GF}(2)$. @@ -271,6 +252,13 @@ class factory :func:`~galois.GF`. galois-fields """ + def astype(self, dtype, **kwargs): + if dtype is GF2BP: + return GF2BP(self) # bits are packed in initialization + + return super().astype(dtype, **kwargs) + + @export class GF2BP( FieldArray, @@ -284,7 +272,7 @@ class GF2BP( bitpacked=True, ): r""" - A :obj:`~galois.FieldArray` subclass over $\mathrm{GF}(2)$. + A :obj:`~galois.FieldArray` subclass over $\mathrm{GF}(2)$ with a bit-packed representation. .. info:: @@ -316,6 +304,29 @@ class factory :func:`~galois.GF`. galois-fields """ + def __init__( + self, x: ElementLike | ArrayLike, dtype: DTypeLike | None = None, **kwargs + ): + if isinstance(x, np.ndarray): + self.view(np.ndarray)[:] = np.packbits(x.view(np.ndarray), axis=-1) + self._unpacked_shape = x.shape + return + + raise RuntimeError("BLAH") + # super().__init__(x, dtype, **kwargs) + + def astype(self, dtype, **kwargs): + if dtype is GF2: + return GF2( + numpy.unpackbits( + self.view(np.ndarray), + axis=-1, + count=self._unpacked_shape[-1], + ) + ) + + return super().astype(dtype, **kwargs) + GF2._default_ufunc_mode = "jit-calculate" GF2._ufunc_modes = ["jit-calculate", "python-calculate"] From 27048bc6bbd893d358a3a1cf05b88a25d3685a62 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Wed, 11 Dec 2024 15:34:59 -0800 Subject: [PATCH 04/21] Check in working matmul (w/ updated numpy) --- requirements-min.txt | 4 +-- src/galois/_fields/_array.py | 26 +++++++++---------- src/galois/_fields/_gf2.py | 48 ++++++++++++++++++++++++------------ 3 files changed, 47 insertions(+), 31 deletions(-) diff --git a/requirements-min.txt b/requirements-min.txt index 92cff6149..2e899ed6a 100644 --- a/requirements-min.txt +++ b/requirements-min.txt @@ -1,13 +1,13 @@ numpy==1.26.0; python_version=="3.12" numpy==1.23.2; python_version=="3.11" numpy==1.21.3; python_version=="3.10" -numpy==1.21.0; python_version=="3.9" +numpy==2.0.0; python_version=="3.9" numpy==1.21.0; python_version=="3.8" numpy==1.21.0; python_version=="3.7" numba==0.59.0; python_version=="3.12" numba==0.57.0; python_version=="3.11" numba==0.55.0; python_version=="3.10" -numba==0.55.0; python_version=="3.9" +numba==0.60.0; python_version=="3.9" numba==0.55.0; python_version=="3.8" numba==0.55.0; python_version=="3.7" typing_extensions==4.0.0 diff --git a/src/galois/_fields/_array.py b/src/galois/_fields/_array.py index c995039dd..9422ba97a 100644 --- a/src/galois/_fields/_array.py +++ b/src/galois/_fields/_array.py @@ -961,19 +961,19 @@ def primitive_roots_of_unity(cls, n: int) -> Self: # Instance methods ############################################################################### - @property - def shape(self) -> npt._shape: - if self._bitpacked and hasattr(self, "_unpacked_shape"): - return self._unpacked_shape - - return super().shape - - @shape.setter - def shape(self, value: npt._shape) -> None: - if self._bitpacked: - raise NotImplementedError("You cannot set the shape on a bit-packed array.") - - super().shape = value + # @property + # def shape(self) -> npt._shape: + # if self._bitpacked and hasattr(self, "_unpacked_shape"): + # return self._unpacked_shape + # + # return super().shape + # + # @shape.setter + # def shape(self, value: npt._shape) -> None: + # if self._bitpacked: + # raise NotImplementedError("You cannot set the shape on a bit-packed array.") + # + # super().shape = value # def __array__(self) -> np.ndarray: # if self._bitpacked and hasattr(self, "_unpacked_shape"): diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 78c49a036..1db95b77e 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -4,10 +4,8 @@ from __future__ import annotations -import numpy import numpy as np -from ..typing import ElementLike, ArrayLike, DTypeLike from .._domains._lookup import ( add_ufunc, divide_ufunc, @@ -21,6 +19,7 @@ ) from .._domains._ufunc import UFuncMixin, matmul_ufunc from .._helper import export +from ..typing import ArrayLike, DTypeLike, ElementLike from ._array import FieldArray @@ -156,19 +155,38 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): a, b = inputs assert isinstance(a, GF2BP) and isinstance(b, GF2BP) - assert a.shape[-1] == b.shape[0] - # a has rows packed, so repack b's columns - b = np.packbits( - np.unpackbits(b.view(np.ndarray), axis=-1, count=b._unpacked_shape[-1]), - axis=0, - ).view(GF2BP) + # a has rows packed, so unpack and repack b's columns + field = self.field + unpacked_shape = b._unpacked_shape + b = field._view( + np.packbits( + np.unpackbits(b.view(np.ndarray), axis=-1, count=b._unpacked_shape[-1]), + axis=0, + ) + ) + b._unpacked_shape = unpacked_shape + + # Make sure the inner dimensions match (e.g. (M, N) x (N, P) -> (M, P)) + assert a.shape[-1] == b.shape[1] + if len(b.shape) == 1: + final_shape = (a.shape[0],) + else: + final_shape = (a.shape[0], b.shape[-1]) inputs = (a, b) - output = super().__call__(ufunc, method, inputs, kwargs, meta) - unpacked_shape = output.shape - output = self.field._view(np.packbits(output.view(np.ndarray), axis=-1)) - output._unpacked_shape = unpacked_shape + # output = super().__call__(ufunc, method, inputs, kwargs, meta) + # output._unpacked_shape = a._unpacked_shape + if len(b.shape) == 1: + output = np.bitwise_xor.reduce(np.unpackbits((a & b).view(np.ndarray), axis=-1), axis=-1) + else: + output = GF2.Zeros(final_shape) + for i in range(b.shape[-1]): + # output[:, i] = np.bitwise_xor.reduce(np.unpackbits((a & b[:, i]).view(np.ndarray), axis=-1), axis=-1) + output[:, i] = np.bitwise_xor.reduce(np.bitwise_count((a & b[:, i]).view(np.ndarray)), axis=-1) % 2 + output = field._view(np.packbits(output.view(np.ndarray), axis=-1)) + # output = field._view(output) + output._unpacked_shape = final_shape return output @@ -304,9 +322,7 @@ class factory :func:`~galois.GF`. galois-fields """ - def __init__( - self, x: ElementLike | ArrayLike, dtype: DTypeLike | None = None, **kwargs - ): + def __init__(self, x: ElementLike | ArrayLike, dtype: DTypeLike | None = None, **kwargs): if isinstance(x, np.ndarray): self.view(np.ndarray)[:] = np.packbits(x.view(np.ndarray), axis=-1) self._unpacked_shape = x.shape @@ -318,7 +334,7 @@ def __init__( def astype(self, dtype, **kwargs): if dtype is GF2: return GF2( - numpy.unpackbits( + np.unpackbits( self.view(np.ndarray), axis=-1, count=self._unpacked_shape[-1], From 3410deb5a041c00e72c89d5aa33ff5d016740dcd Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Wed, 11 Dec 2024 15:39:12 -0800 Subject: [PATCH 05/21] Clean up matmul function --- src/galois/_fields/_gf2.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 1db95b77e..f7ec103c6 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -156,7 +156,7 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): assert isinstance(a, GF2BP) and isinstance(b, GF2BP) - # a has rows packed, so unpack and repack b's columns + # bit-packed matrices have rows packed by default, so unpack the second operand and repack to columns field = self.field unpacked_shape = b._unpacked_shape b = field._view( @@ -174,18 +174,15 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): else: final_shape = (a.shape[0], b.shape[-1]) - inputs = (a, b) - # output = super().__call__(ufunc, method, inputs, kwargs, meta) - # output._unpacked_shape = a._unpacked_shape if len(b.shape) == 1: + # matrix-vector multiplication output = np.bitwise_xor.reduce(np.unpackbits((a & b).view(np.ndarray), axis=-1), axis=-1) else: + # matrix-matrix multiplication output = GF2.Zeros(final_shape) for i in range(b.shape[-1]): - # output[:, i] = np.bitwise_xor.reduce(np.unpackbits((a & b[:, i]).view(np.ndarray), axis=-1), axis=-1) output[:, i] = np.bitwise_xor.reduce(np.bitwise_count((a & b[:, i]).view(np.ndarray)), axis=-1) % 2 output = field._view(np.packbits(output.view(np.ndarray), axis=-1)) - # output = field._view(output) output._unpacked_shape = final_shape return output From d9ba03cab5cc2f7598bbd64c7cef01486f032a71 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Wed, 11 Dec 2024 15:43:00 -0800 Subject: [PATCH 06/21] Remove class-level bitpacked field --- src/galois/_domains/_array.py | 73 +++----------------------------- src/galois/_domains/_function.py | 6 +-- src/galois/_domains/_meta.py | 3 -- src/galois/_fields/_array.py | 42 ------------------ src/galois/_fields/_gf2.py | 1 - 5 files changed, 8 insertions(+), 117 deletions(-) diff --git a/src/galois/_domains/_array.py b/src/galois/_domains/_array.py index 721eacb35..283de9d49 100644 --- a/src/galois/_domains/_array.py +++ b/src/galois/_domains/_array.py @@ -7,7 +7,6 @@ import abc import contextlib import random -from numbers import Integral from typing import Generator import numpy as np @@ -54,19 +53,10 @@ def __new__( dtype = cls._get_dtype(dtype) x = cls._verify_array_like_types_and_values(x) - packed = False - if cls._bitpacked and isinstance(x, np.ndarray) and x.size > 0: - array = np.packbits(x.view(np.ndarray), axis=-1) - packed = True - else: - array = np.array(x, dtype=dtype, copy=copy, order=order, ndmin=ndmin) + array = np.array(x, dtype=dtype, copy=copy, order=order, ndmin=ndmin) # Perform view without verification since the elements were verified in _verify_array_like_types_and_values() - galois_array = cls._view(array) - if cls._bitpacked and packed: - galois_array._unpacked_shape = x.shape - galois_array.packed = packed - return galois_array + return cls._view(array) @classmethod def _get_dtype(cls, dtype: DTypeLike) -> np.dtype: @@ -182,17 +172,7 @@ def Zeros(cls, shape: ShapeLike, dtype: DTypeLike | None = None) -> Self: """ dtype = cls._get_dtype(dtype) array = np.zeros(shape, dtype=dtype) - if cls._bitpacked: - original_shape = array.shape - array = np.packbits(array, axis=-1) - - array = cls._view(array) - - if cls._bitpacked: - array.packed = True - array._unpacked_shape = original_shape - - return array + return cls._view(array) @classmethod def Ones(cls, shape: ShapeLike, dtype: DTypeLike | None = None) -> Self: @@ -210,17 +190,7 @@ def Ones(cls, shape: ShapeLike, dtype: DTypeLike | None = None) -> Self: """ dtype = cls._get_dtype(dtype) array = np.ones(shape, dtype=dtype) - if cls._bitpacked: - original_shape = array.shape - array = np.packbits(array, axis=-1) - - array = cls._view(array) - - if cls._bitpacked: - array.packed = True - array._unpacked_shape = original_shape - - return array + return cls._view(array) @classmethod def Range( @@ -257,17 +227,8 @@ def Range( raise ValueError(f"Argument 'stop' must be within the field's order {cls.order}, not {stop}.") array = np.arange(start, stop, step=step, dtype=dtype) - if cls._bitpacked: - original_shape = array.shape - array = np.packbits(array, axis=-1) - - array = cls._view(array) - if cls._bitpacked: - array.packed = True - array._unpacked_shape = original_shape - - return array + return cls._view(array) @classmethod def Random( @@ -334,17 +295,7 @@ def Random( for _ in iterator: array[iterator.multi_index] = random.randint(low, high - 1) - if cls._bitpacked: - original_shape = array.shape - array = np.packbits(array, axis=-1) - - array = cls._view(array) - - if cls._bitpacked: - array.packed = True - array._unpacked_shape = original_shape - - return array + return cls._view(array) @classmethod def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: @@ -362,17 +313,7 @@ def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: """ dtype = cls._get_dtype(dtype) array = np.identity(size, dtype=dtype) - if cls._bitpacked: - original_shape = array.shape - array = np.packbits(array, axis=-1) - - array = cls._view(array) - - if cls._bitpacked: - array.packed = True - array._unpacked_shape = original_shape - - return array + return cls._view(array) ############################################################################### # Ufunc compilation routines diff --git a/src/galois/_domains/_function.py b/src/galois/_domains/_function.py index 30c98d576..eaf7d4a9f 100644 --- a/src/galois/_domains/_function.py +++ b/src/galois/_domains/_function.py @@ -334,7 +334,6 @@ class FunctionMixin(np.ndarray, metaclass=ArrayMeta): np.convolve: "_convolve", np.fft.fft: "_fft", np.fft.ifft: "_ifft", - np.array_equal: "_array_equal", } _convolve: Function @@ -354,10 +353,7 @@ def __array_function__(self, func, types, args, kwargs): field = type(self) if func in field._OVERRIDDEN_FUNCTIONS: - try: - output = getattr(field, field._OVERRIDDEN_FUNCTIONS[func])(*args, **kwargs) - except AttributeError: - output = super().__array_function__(func, types, args, kwargs) + output = getattr(field, field._OVERRIDDEN_FUNCTIONS[func])(*args, **kwargs) elif func in field._UNSUPPORTED_FUNCTIONS: raise NotImplementedError( diff --git a/src/galois/_domains/_meta.py b/src/galois/_domains/_meta.py index 85d726975..892bf549f 100644 --- a/src/galois/_domains/_meta.py +++ b/src/galois/_domains/_meta.py @@ -35,9 +35,6 @@ def __init__(cls, name, bases, namespace, **kwargs): cls._irreducible_poly_int: int = kwargs.get("irreducible_poly_int", 0) cls._primitive_element: int = kwargs.get("primitive_element", 0) cls._dtypes = cls._determine_dtypes() - cls._bitpacked = kwargs.get("bitpacked", False) - if cls._bitpacked and (cls._order != 2 or cls._degree != 1): - raise ValueError("Only fields over GF(2) support bit-packing.") if cls._dtypes == [np.object_]: cls._default_ufunc_mode = "python-calculate" diff --git a/src/galois/_fields/_array.py b/src/galois/_fields/_array.py index 9422ba97a..3a074bc92 100644 --- a/src/galois/_fields/_array.py +++ b/src/galois/_fields/_array.py @@ -7,7 +7,6 @@ from typing import Generator import numpy as np -import numpy.typing as npt from typing_extensions import Literal, Self from .._domains import Array, _linalg @@ -961,47 +960,6 @@ def primitive_roots_of_unity(cls, n: int) -> Self: # Instance methods ############################################################################### - # @property - # def shape(self) -> npt._shape: - # if self._bitpacked and hasattr(self, "_unpacked_shape"): - # return self._unpacked_shape - # - # return super().shape - # - # @shape.setter - # def shape(self, value: npt._shape) -> None: - # if self._bitpacked: - # raise NotImplementedError("You cannot set the shape on a bit-packed array.") - # - # super().shape = value - - # def __array__(self) -> np.ndarray: - # if self._bitpacked and hasattr(self, "_unpacked_shape"): - # return np.unpackbits(self.view(np.ndarray), axis=-1, count=self._unpacked_shape[-1]) - # - # return self - # - # def __getitem__(self, index): - # return self.__array__()[index] - # - # def __setitem__(self, index, value): - # array_unpacked = self.__array__() - # array_unpacked[index] = value - # self.view(np.ndarray)[:] = np.packbits(array_unpacked, axis=-1)[:] - - @property - def T(self) -> Self: - if self._bitpacked and hasattr(self, "_unpacked_shape"): - transposed = np.unpackbits( - self.view(np.ndarray), axis=-1, count=self._unpacked_shape[-1] - ).T - original_shape = transposed.shape - transposed = np.packbits(transposed, axis=-1) - transposed.original_shape = original_shape - return transposed - - return super().T - def additive_order(self) -> int | np.ndarray: r""" Computes the additive order of each element in $x$. diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index f7ec103c6..03c2da573 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -284,7 +284,6 @@ class GF2BP( irreducible_poly_int=3, is_primitive_poly=True, primitive_element=1, - bitpacked=True, ): r""" A :obj:`~galois.FieldArray` subclass over $\mathrm{GF}(2)$ with a bit-packed representation. From 7bf44e4d24e3b3675a0f178e2cea45785d40543c Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Wed, 11 Dec 2024 17:40:16 -0800 Subject: [PATCH 07/21] Fix matmul shape bug and __new__ --- src/galois/_fields/_gf2.py | 30 +++++++++++++++++++++++------- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 03c2da573..2ef21c5c3 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -5,6 +5,7 @@ from __future__ import annotations import numpy as np +from typing_extensions import Literal, Self from .._domains._lookup import ( add_ufunc, @@ -168,7 +169,7 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): b._unpacked_shape = unpacked_shape # Make sure the inner dimensions match (e.g. (M, N) x (N, P) -> (M, P)) - assert a.shape[-1] == b.shape[1] + assert a.shape[-1] == b.shape[0] if len(b.shape) == 1: final_shape = (a.shape[0],) else: @@ -318,14 +319,29 @@ class factory :func:`~galois.GF`. galois-fields """ - def __init__(self, x: ElementLike | ArrayLike, dtype: DTypeLike | None = None, **kwargs): + def __new__( + cls, + x: ElementLike | ArrayLike, + dtype: DTypeLike | None = None, + copy: bool = True, + order: Literal["K", "A", "C", "F"] = "K", + ndmin: int = 0, + ) -> Self: if isinstance(x, np.ndarray): - self.view(np.ndarray)[:] = np.packbits(x.view(np.ndarray), axis=-1) - self._unpacked_shape = x.shape - return + dtype = cls._get_dtype(dtype) - raise RuntimeError("BLAH") - # super().__init__(x, dtype, **kwargs) + x = cls._verify_array_like_types_and_values(x) + array = cls._view(np.packbits(np.array(x, dtype=dtype, copy=copy, order=order, ndmin=ndmin).view(np.ndarray), axis=-1)) + array._unpacked_shape = x.shape + + # Perform view without verification since the elements were verified in _verify_array_like_types_and_values() + return array + + raise NotImplementedError( + "GF2BP is a custom bit-packed GF2 class with limited functionality. " + "If you were using an alternate constructor (e.g. Random), then use the GF2 class and convert it to the " + "bit-packed version by using `.astype(GF2BP)`." + ) def astype(self, dtype, **kwargs): if dtype is GF2: From 7054295b90e7bda9737916a7ea069928113534e8 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Fri, 13 Dec 2024 16:52:31 -0800 Subject: [PATCH 08/21] Override np.packbits and np.unpackbits --- src/galois/_domains/_function.py | 9 ++- src/galois/_fields/_gf2.py | 96 ++++++++++++++++++++------------ 2 files changed, 66 insertions(+), 39 deletions(-) diff --git a/src/galois/_domains/_function.py b/src/galois/_domains/_function.py index eaf7d4a9f..c206cf707 100644 --- a/src/galois/_domains/_function.py +++ b/src/galois/_domains/_function.py @@ -304,8 +304,6 @@ class FunctionMixin(np.ndarray, metaclass=ArrayMeta): _UNSUPPORTED_FUNCTIONS = [ # Unary - np.packbits, - np.unpackbits, np.unwrap, np.around, np.round, @@ -334,6 +332,8 @@ class FunctionMixin(np.ndarray, metaclass=ArrayMeta): np.convolve: "_convolve", np.fft.fft: "_fft", np.fft.ifft: "_ifft", + np.packbits: "_packbits", + np.unpackbits: "_unpackbits", } _convolve: Function @@ -353,7 +353,10 @@ def __array_function__(self, func, types, args, kwargs): field = type(self) if func in field._OVERRIDDEN_FUNCTIONS: - output = getattr(field, field._OVERRIDDEN_FUNCTIONS[func])(*args, **kwargs) + try: + output = getattr(field, field._OVERRIDDEN_FUNCTIONS[func])(*args, **kwargs) + except AttributeError: + output = super().__array_function__(func, types, args, kwargs) elif func in field._UNSUPPORTED_FUNCTIONS: raise NotImplementedError( diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 2ef21c5c3..72e91659e 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -5,7 +5,7 @@ from __future__ import annotations import numpy as np -from typing_extensions import Literal, Self +from typing_extensions import Literal, Self, Optional from .._domains._lookup import ( add_ufunc, @@ -84,6 +84,30 @@ class sqrt(sqrt_ufunc): def implementation(self, a: FieldArray) -> FieldArray: return a.copy() +def packbits(a, axis=None, bitorder='big'): + if isinstance(a, GF2BP): + return a + + if not isinstance(a, GF2): + raise TypeError("Bit-packing is only supported on instances of GF2.") + + axis = -1 if axis is None else axis + packed = GF2BP(np.packbits(a.view(np.ndarray), axis=axis, bitorder=bitorder), a.shape[axis]) + return packed + + +def unpackbits(a, axis=None, count=None, bitorder='big'): + if isinstance(a, GF2): + return a + + if not isinstance(a, GF2BP): + raise TypeError("Bit-unpacking is only supported on instances of GF2BP.") + + if axis is None: + axis = -1 + + return GF2(np.unpackbits(a.view(np.ndarray), axis=axis, count=a._axis_count if count is None else count, bitorder=bitorder)) + class UFuncMixin_2_1(UFuncMixin): """ @@ -101,6 +125,8 @@ def __init_subclass__(cls) -> None: cls._power = power(cls) cls._log = log(cls) cls._sqrt = sqrt(cls) + cls._packbits = packbits + cls._unpackbits = unpackbits class add_ufunc_bitpacked(add_ufunc): @@ -110,7 +136,7 @@ class add_ufunc_bitpacked(add_ufunc): def __call__(self, ufunc, method, inputs, kwargs, meta): output = super().__call__(ufunc, method, inputs, kwargs, meta) - output._unpacked_shape = inputs[0]._unpacked_shape + output._axis_count = inputs[0]._axis_count return output @@ -121,7 +147,7 @@ class subtract_ufunc_bitpacked(subtract_ufunc): def __call__(self, ufunc, method, inputs, kwargs, meta): output = super().__call__(ufunc, method, inputs, kwargs, meta) - output._unpacked_shape = inputs[0]._unpacked_shape + output._axis_count = inputs[0]._axis_count return output @@ -132,7 +158,7 @@ class multiply_ufunc_bitpacked(multiply_ufunc): def __call__(self, ufunc, method, inputs, kwargs, meta): output = super().__call__(ufunc, method, inputs, kwargs, meta) - output._unpacked_shape = inputs[0]._unpacked_shape + output._axis_count = inputs[0]._axis_count return output @@ -143,7 +169,7 @@ class divide_ufunc_bitpacked(divide): def __call__(self, ufunc, method, inputs, kwargs, meta): output = super().__call__(ufunc, method, inputs, kwargs, meta) - output._unpacked_shape = inputs[0]._unpacked_shape + output._axis_count = inputs[0]._axis_count return output @@ -157,16 +183,16 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): assert isinstance(a, GF2BP) and isinstance(b, GF2BP) - # bit-packed matrices have rows packed by default, so unpack the second operand and repack to columns + # bit-packed matrices have columns packed by default, so unpack the second operand and repack to rows field = self.field - unpacked_shape = b._unpacked_shape + row_axis_count = b.shape[0] b = field._view( np.packbits( - np.unpackbits(b.view(np.ndarray), axis=-1, count=b._unpacked_shape[-1]), + np.unpackbits(b.view(np.ndarray), axis=-1, count=b._axis_count), axis=0, ) ) - b._unpacked_shape = unpacked_shape + b._axis_count = row_axis_count # Make sure the inner dimensions match (e.g. (M, N) x (N, P) -> (M, P)) assert a.shape[-1] == b.shape[0] @@ -182,9 +208,11 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): # matrix-matrix multiplication output = GF2.Zeros(final_shape) for i in range(b.shape[-1]): + # TODO: Include alternate path for numpy < v2 + # output[:, i] = np.bitwise_xor.reduce(np.unpackbits((a & b[:, i]).view(np.ndarray), axis=-1), axis=-1) output[:, i] = np.bitwise_xor.reduce(np.bitwise_count((a & b[:, i]).view(np.ndarray)), axis=-1) % 2 output = field._view(np.packbits(output.view(np.ndarray), axis=-1)) - output._unpacked_shape = final_shape + output._axis_count = final_shape[-1] return output @@ -192,7 +220,6 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): def not_implemented(*args, **kwargs): return NotImplemented - class UFuncMixin_2_1_BitPacked(UFuncMixin): """ A mixin class that provides explicit calculation arithmetic for GF(2). @@ -209,6 +236,8 @@ def __init_subclass__(cls) -> None: cls._power = power(cls) cls._log = log(cls) cls._sqrt = sqrt(cls) + cls._packbits = packbits + cls._unpackbits = unpackbits @classmethod def _assign_ufuncs(cls): @@ -268,12 +297,6 @@ class factory :func:`~galois.GF`. galois-fields """ - def astype(self, dtype, **kwargs): - if dtype is GF2BP: - return GF2BP(self) # bits are packed in initialization - - return super().astype(dtype, **kwargs) - @export class GF2BP( @@ -322,39 +345,40 @@ class factory :func:`~galois.GF`. def __new__( cls, x: ElementLike | ArrayLike, + axis_element_count: Optional[int] = None, dtype: DTypeLike | None = None, copy: bool = True, order: Literal["K", "A", "C", "F"] = "K", ndmin: int = 0, ) -> Self: - if isinstance(x, np.ndarray): - dtype = cls._get_dtype(dtype) + # axis_element_count is required, but by making it optional it allows us to catch uses of the class that are not + # supported (e.g. Random) + if isinstance(x, np.ndarray) and axis_element_count is not None: + # NOTE: I'm not sure that we want to change the dtype specifically for the bit-packed version or how we verify + # dtype = cls._get_dtype(dtype) + # x = cls._verify_array_like_types_and_values(x) - x = cls._verify_array_like_types_and_values(x) - array = cls._view(np.packbits(np.array(x, dtype=dtype, copy=copy, order=order, ndmin=ndmin).view(np.ndarray), axis=-1)) - array._unpacked_shape = x.shape + array = cls._view(np.array(x, dtype=dtype, copy=copy, order=order, ndmin=ndmin)) + array._axis_count = axis_element_count - # Perform view without verification since the elements were verified in _verify_array_like_types_and_values() return array raise NotImplementedError( "GF2BP is a custom bit-packed GF2 class with limited functionality. " "If you were using an alternate constructor (e.g. Random), then use the GF2 class and convert it to the " - "bit-packed version by using `.astype(GF2BP)`." + "bit-packed version by using `np.packbits`." ) - def astype(self, dtype, **kwargs): - if dtype is GF2: - return GF2( - np.unpackbits( - self.view(np.ndarray), - axis=-1, - count=self._unpacked_shape[-1], - ) - ) - - return super().astype(dtype, **kwargs) - + def __init__( + self, + x: ElementLike | ArrayLike, + axis_element_count: Optional[int] = None, + dtype: DTypeLike | None = None, + copy: bool = True, + order: Literal["K", "A", "C", "F"] = "K", + ndmin: int = 0, + ): + pass GF2._default_ufunc_mode = "jit-calculate" GF2._ufunc_modes = ["jit-calculate", "python-calculate"] From b9d75e02e172147abfa0bcf4b5f94e0956b23c43 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Fri, 13 Dec 2024 17:33:09 -0800 Subject: [PATCH 09/21] Disallow unsupported numpy functions --- src/galois/_fields/_gf2.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 72e91659e..785151b41 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -147,7 +147,7 @@ class subtract_ufunc_bitpacked(subtract_ufunc): def __call__(self, ufunc, method, inputs, kwargs, meta): output = super().__call__(ufunc, method, inputs, kwargs, meta) - output._axis_count = inputs[0]._axis_count + output._axis_count = max(i._axis_count for i in inputs) return output @@ -158,7 +158,7 @@ class multiply_ufunc_bitpacked(multiply_ufunc): def __call__(self, ufunc, method, inputs, kwargs, meta): output = super().__call__(ufunc, method, inputs, kwargs, meta) - output._axis_count = inputs[0]._axis_count + output._axis_count = max(i._axis_count for i in inputs) return output @@ -169,7 +169,7 @@ class divide_ufunc_bitpacked(divide): def __call__(self, ufunc, method, inputs, kwargs, meta): output = super().__call__(ufunc, method, inputs, kwargs, meta) - output._axis_count = inputs[0]._axis_count + output._axis_count = max(i._axis_count for i in inputs) return output @@ -218,6 +218,7 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): def not_implemented(*args, **kwargs): + # TODO: Add a better error message about limited support return NotImplemented class UFuncMixin_2_1_BitPacked(UFuncMixin): @@ -231,10 +232,10 @@ def __init_subclass__(cls) -> None: cls._negative = negative_ufunc(cls, override=np.positive) cls._subtract = subtract_ufunc_bitpacked(cls, override=np.bitwise_xor) cls._multiply = multiply_ufunc_bitpacked(cls, override=np.bitwise_and) - cls._reciprocal = reciprocal(cls) + cls._reciprocal = not_implemented # reciprocal(cls) cls._divide = divide_ufunc_bitpacked(cls) - cls._power = power(cls) - cls._log = log(cls) + cls._power = not_implemented # power(cls) + cls._log = not_implemented # log(cls) cls._sqrt = sqrt(cls) cls._packbits = packbits cls._unpackbits = unpackbits From 1186fa4cb5eaf6286c9f1158b51a1ab277c97602 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Thu, 19 Dec 2024 16:42:29 -0800 Subject: [PATCH 10/21] Start work on concatenate, inv, and indexing --- src/galois/_domains/_function.py | 5 +- src/galois/_fields/_gf2.py | 106 +++++++++++++++++++++++++++++-- 2 files changed, 105 insertions(+), 6 deletions(-) diff --git a/src/galois/_domains/_function.py b/src/galois/_domains/_function.py index c206cf707..f7f809b5f 100644 --- a/src/galois/_domains/_function.py +++ b/src/galois/_domains/_function.py @@ -334,6 +334,7 @@ class FunctionMixin(np.ndarray, metaclass=ArrayMeta): np.fft.ifft: "_ifft", np.packbits: "_packbits", np.unpackbits: "_unpackbits", + np.concatenate: "_concatenate", } _convolve: Function @@ -375,8 +376,8 @@ def __array_function__(self, func, types, args, kwargs): output = super().__array_function__(func, types, args, kwargs) - if func in field._FUNCTIONS_REQUIRING_VIEW: - output = field._view(output) if not np.isscalar(output) else field(output, dtype=self.dtype) + if func in field._FUNCTIONS_REQUIRING_VIEW: + output = field._view(output) if not np.isscalar(output) else field(output, dtype=self.dtype) return output diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 785151b41..b3d926647 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -4,9 +4,12 @@ from __future__ import annotations +from typing import Sequence + import numpy as np from typing_extensions import Literal, Self, Optional +from .._domains._linalg import row_reduce_jit from .._domains._lookup import ( add_ufunc, divide_ufunc, @@ -19,7 +22,9 @@ subtract_ufunc, ) from .._domains._ufunc import UFuncMixin, matmul_ufunc -from .._helper import export +from .._domains._function import Function +from .._domains._array import Array +from .._helper import export, verify_isinstance from ..typing import ArrayLike, DTypeLike, ElementLike from ._array import FieldArray @@ -216,6 +221,50 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): return output +class concatenate_bitpacked(Function): + """Concatenates matrices together""" + + def __call__(self, arrays: list[Array], **kwargs): + unpacked_arrays = [] + for array in arrays: + verify_isinstance(array, self.field) + unpacked_arrays.append(np.unpackbits(array)) + + unpacked = np.concatenate(unpacked_arrays, **kwargs) + return np.packbits(unpacked) + + +class inv_bitpacked(Function): + """ + Computes the inverse of the square matrix. + """ + + def __call__(self, A: Array) -> Array: + verify_isinstance(A, self.field) + # if not (A.ndim == 2 and A.shape[0] == A.shape[1]): + # raise np.linalg.LinAlgError(f"Argument 'A' must be square, not {A.shape}.") + + n = A.shape[0] + I = self.field.Identity(n, dtype=A.dtype) + + # Concatenate A and I to get the matrix AI = [A | I] + AI = np.concatenate((A, I), axis=-1) + + # Perform Gaussian elimination to get the reduced row echelon form AI_rre = [I | A^-1] + AI_rre, _ = row_reduce_jit(self.field)(AI, ncols=n) + + # The rank is the number of non-zero rows of the row reduced echelon form + rank = np.sum(~np.all(AI_rre[:, 0:n] == 0, axis=1)) + if not rank == n: + raise np.linalg.LinAlgError( + f"Argument 'A' is singular and not invertible because it does not have full rank of {n}, " + f"but rank of {rank}." + ) + + A_inv = AI_rre[:, -n:] + + return A_inv + def not_implemented(*args, **kwargs): # TODO: Add a better error message about limited support @@ -239,6 +288,7 @@ def __init_subclass__(cls) -> None: cls._sqrt = sqrt(cls) cls._packbits = packbits cls._unpackbits = unpackbits + cls._concatenate = concatenate_bitpacked(cls) @classmethod def _assign_ufuncs(cls): @@ -246,6 +296,7 @@ def _assign_ufuncs(cls): # We have to set this here because ArrayMeta would override it. cls._matmul = matmul_ufunc_bitpacked(cls) + cls._inv = inv_bitpacked(cls) # NOTE: There is a "verbatim" block in the docstring because we were not able to monkey-patch GF2 like the @@ -381,10 +432,57 @@ def __init__( ): pass + @classmethod + def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: + r""" + Creates an $n \times n$ identity matrix. + + Arguments: + size: The size $n$ along one dimension of the identity matrix. + dtype: The :obj:`numpy.dtype` of the array elements. The default is `None` which represents the smallest + unsigned data type for this :obj:`~galois.Array` subclass (the first element in + :obj:`~galois.Array.dtypes`). + + Returns: + A 2-D identity matrix with shape `(size, size)`. + """ + array = GF2.Identity(size, dtype=dtype) + return np.packbits(array) + + def get_unpacked_slice(self, index): + if isinstance(index, Sequence): + if len(index) == 2: + row_index, col_index = index + if isinstance(col_index, int): + post_index = (slice(None), col_index) + index = (row_index, col_index // 8) + elif isinstance(col_index, slice): + post_index = (slice(None), col_index) + col_index = slice(col_index.start // 8, max(col_index.step // 8, 1), max(col_index.stop // 8, 1)) + index = (row_index, col_index) + elif isinstance(col_index, Sequence): + post_index = (list(range(len(row_index))), col_index) + col_index = tuple(s // 8 for s in col_index) + index = (row_index, col_index) + elif isinstance(index, slice): + # TODO + pass + + packed = self[index] + if len(packed.shape) == 1: + packed = packed[:, None] + unpacked = np.unpackbits(packed, axis=-1, count=self._axis_count) + return GF2._view(unpacked[post_index]) + + + def set_unpacked_slice(self, index, value): + pass + + GF2._default_ufunc_mode = "jit-calculate" GF2._ufunc_modes = ["jit-calculate", "python-calculate"] GF2.compile("auto") -GF2BP._default_ufunc_mode = "jit-calculate" -GF2BP._ufunc_modes = ["jit-calculate", "python-calculate"] -GF2BP.compile("auto") +# GF2BP._default_ufunc_mode = "jit-calculate" +# GF2BP._ufunc_modes = ["jit-calculate", "python-calculate"] +# GF2BP.compile("auto") From c07ba0f462d5171da3256e3a698953c5deb5b0b5 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Fri, 20 Dec 2024 16:59:32 -0800 Subject: [PATCH 11/21] Add additional indexing support --- src/galois/_fields/_gf2.py | 48 +++++++++++++++++++++---- tests/fields/test_bitpacked.py | 65 ++++++++++++++++++++++++++++++++++ 2 files changed, 106 insertions(+), 7 deletions(-) create mode 100644 tests/fields/test_bitpacked.py diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index b3d926647..027e88e2d 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -405,7 +405,7 @@ def __new__( ) -> Self: # axis_element_count is required, but by making it optional it allows us to catch uses of the class that are not # supported (e.g. Random) - if isinstance(x, np.ndarray) and axis_element_count is not None: + if isinstance(x, (tuple, list, np.ndarray, FieldArray)) and axis_element_count is not None: # NOTE: I'm not sure that we want to change the dtype specifically for the bit-packed version or how we verify # dtype = cls._get_dtype(dtype) # x = cls._verify_array_like_types_and_values(x) @@ -450,7 +450,8 @@ def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: return np.packbits(array) def get_unpacked_slice(self, index): - if isinstance(index, Sequence): + post_index = NotImplemented + if isinstance(index, (Sequence, np.ndarray)): if len(index) == 2: row_index, col_index = index if isinstance(col_index, int): @@ -464,16 +465,49 @@ def get_unpacked_slice(self, index): post_index = (list(range(len(row_index))), col_index) col_index = tuple(s // 8 for s in col_index) index = (row_index, col_index) + elif col_index is None: # new axis + post_index = (slice(None), None) + index = (row_index,) + elif ((isinstance(index, np.ndarray) and index.ndim == 1) or + (isinstance(index, list) and all(isinstance(x, int) for x in index))): + post_index = index + index = list(range((len(index) // 8) + 1)) + elif isinstance(index, tuple) and any(x is Ellipsis for x in index): + post_index = index[1:] + axis_adjustment = (slice(None),) if index[-1] is Ellipsis else (index[-1] // 8,) + index = index[:-1] + axis_adjustment elif isinstance(index, slice): - # TODO - pass - - packed = self[index] - if len(packed.shape) == 1: + if self.ndim > 1: + # Rows aren't packed, so we can index normally + post_index = slice(None) + if len(self.shape) == 1: + # Array is 1-D, so we need to adjust + post_index = index + index = slice(index.start // 8 if index.start is not None else index.start, + max(index.step // 8, 1) if index.step is not None else index.step, + max(index.stop // 8, 1) if index.stop is not None else index.stop) + elif isinstance(index, int): + post_index = index + index //= 8 + + if post_index is NotImplemented: + raise NotImplementedError(f"The following indexing scheme is not supported:\n{index}\n" + "If you believe this scheme should be supported, " + "please submit a GitHub issue at https://github.com/mhostetter/galois/issues.\n\n" + "If you'd like to perform this operation on the data, you should first call " + "`array = array.view(np.ndarray)` and then call the function." + ) + + packed = self.view(np.ndarray)[index] + if np.isscalar(packed): + packed = GF2BP([packed], self._axis_count).view(np.ndarray) + if packed.ndim == 1 and self.ndim > 1: packed = packed[:, None] unpacked = np.unpackbits(packed, axis=-1, count=self._axis_count) return GF2._view(unpacked[post_index]) + def __getitem__(self, item): + return self.get_unpacked_slice(item) def set_unpacked_slice(self, index, value): pass diff --git a/tests/fields/test_bitpacked.py b/tests/fields/test_bitpacked.py new file mode 100644 index 000000000..b11cc2336 --- /dev/null +++ b/tests/fields/test_bitpacked.py @@ -0,0 +1,65 @@ +import numpy as np +import galois + +def test_galois_array_indexing(): + # Define a Galois field array + GF = galois.GF(2) + arr = GF([1, 0, 1, 1]) + arr = np.packbits(arr) + + # 1. Basic Indexing + assert arr[0] == GF(1) + assert arr[2] == GF(1) + + # 2. Negative Indexing + assert arr[-1] == GF(1) + assert arr[-2] == GF(1) + + # 3. Slicing + assert np.array_equal(arr[1:3], GF([0, 1])) + assert np.array_equal(arr[:3], GF([1, 0, 1])) + assert np.array_equal(arr[::2], GF([1, 1])) + assert np.array_equal(arr[::-1], GF([1, 1, 0, 1])) + + # 4. Multidimensional Indexing + arr_2d = GF([[1, 0], [0, 1]]) + arr_2d = np.packbits(arr_2d) + assert arr_2d[0, 1] == GF(0) + assert np.array_equal(arr_2d[:, 1], GF([0, 1])) + + # 5. Boolean Indexing + mask = np.array([True, False, True, False]) + assert np.array_equal(arr[mask], GF([1, 1])) + + # 6. Fancy Indexing + indices = [0, 2, 3] + assert np.array_equal(arr[indices], GF([1, 1, 1])) + + # 7. Ellipsis + arr_3d = GF(np.random.randint(0, 2, (2, 3, 4))) + arr_3d = np.packbits(arr_3d) + shape_check = arr_3d[0, ..., 1].shape # (3,) + assert shape_check == (3,) + + # 8. Indexing with slice objects + s = slice(1, 3) + assert np.array_equal(arr[s], GF([0, 1])) + + # 9. Using np.newaxis + reshaped = arr[:, np.newaxis] + assert reshaped.shape == (4, 1) + + # 10. Indexing with np.ix_ + row_indices = np.array([0, 1]) + col_indices = np.array([0, 1]) + sub_matrix = arr_2d[np.ix_(row_indices, col_indices)] + assert np.array_equal(sub_matrix, GF([[1, 0], [0, 1]])) + + # 11. Indexing with np.take + taken = np.take(arr, [0, 2]) + assert np.array_equal(taken, GF([1, 1])) + + print("All tests passed.") + +if __name__ == "__main__": + test_galois_array_indexing() From 7c15a8503e23aad1661b3cc1eae8f36bc45c3af1 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Tue, 24 Dec 2024 07:31:32 -0800 Subject: [PATCH 12/21] Add tests for setitem and implement --- src/galois/_domains/_function.py | 14 +++--- src/galois/_fields/_gf2.py | 43 +++++++++++++++-- tests/fields/test_bitpacked.py | 81 +++++++++++++++++++++++++++++++- 3 files changed, 125 insertions(+), 13 deletions(-) diff --git a/src/galois/_domains/_function.py b/src/galois/_domains/_function.py index f7f809b5f..d551ddbb9 100644 --- a/src/galois/_domains/_function.py +++ b/src/galois/_domains/_function.py @@ -352,12 +352,14 @@ def __array_function__(self, func, types, args, kwargs): Override the standard NumPy function calls with the new finite field functions. """ field = type(self) + output = None if func in field._OVERRIDDEN_FUNCTIONS: try: output = getattr(field, field._OVERRIDDEN_FUNCTIONS[func])(*args, **kwargs) except AttributeError: - output = super().__array_function__(func, types, args, kwargs) + # fall through to use the default numpy function + pass elif func in field._UNSUPPORTED_FUNCTIONS: raise NotImplementedError( @@ -368,12 +370,12 @@ def __array_function__(self, func, types, args, kwargs): "`array = array.view(np.ndarray)` and then call the function." ) - else: - if func is np.insert: - args = list(args) - args[2] = self._verify_array_like_types_and_values(args[2]) - args = tuple(args) + if func is np.insert: + args = list(args) + args[2] = self._verify_array_like_types_and_values(args[2]) + args = tuple(args) + if output is None: output = super().__array_function__(func, types, args, kwargs) if func in field._FUNCTIONS_REQUIRING_VIEW: diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 027e88e2d..0a9c8ec05 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -97,7 +97,8 @@ def packbits(a, axis=None, bitorder='big'): raise TypeError("Bit-packing is only supported on instances of GF2.") axis = -1 if axis is None else axis - packed = GF2BP(np.packbits(a.view(np.ndarray), axis=axis, bitorder=bitorder), a.shape[axis]) + axis_element_count = 1 if a.ndim == 0 else a.shape[axis] + packed = GF2BP(np.packbits(a.view(np.ndarray), axis=axis, bitorder=bitorder), axis_element_count) return packed @@ -449,7 +450,7 @@ def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: array = GF2.Identity(size, dtype=dtype) return np.packbits(array) - def get_unpacked_slice(self, index): + def get_index_parameters(self, index): post_index = NotImplemented if isinstance(index, (Sequence, np.ndarray)): if len(index) == 2: @@ -461,8 +462,11 @@ def get_unpacked_slice(self, index): post_index = (slice(None), col_index) col_index = slice(col_index.start // 8, max(col_index.step // 8, 1), max(col_index.stop // 8, 1)) index = (row_index, col_index) - elif isinstance(col_index, Sequence): - post_index = (list(range(len(row_index))), col_index) + elif isinstance(col_index, (Sequence, np.ndarray)): + if isinstance(row_index, np.ndarray): + post_index = np.array(range(len(row_index))).reshape(row_index.shape), col_index + else: + post_index = list(range(len(row_index))), col_index col_index = tuple(s // 8 for s in col_index) index = (row_index, col_index) elif col_index is None: # new axis @@ -476,6 +480,13 @@ def get_unpacked_slice(self, index): post_index = index[1:] axis_adjustment = (slice(None),) if index[-1] is Ellipsis else (index[-1] // 8,) index = index[:-1] + axis_adjustment + elif isinstance(index, tuple) and any(isinstance(x, slice) for x in index): + post_index = index[1:] + axis_adjustment = (slice(index.start // 8 if index.start is not None else index.start, + max(index.step // 8, 1) if index.step is not None else index.step, + max(index.stop // 8, 1) if index.stop is not None else index.stop) + if isinstance(index[-1], slice) else (index[-1] // 8,)) + index = index[:-1] + axis_adjustment elif isinstance(index, slice): if self.ndim > 1: # Rows aren't packed, so we can index normally @@ -490,6 +501,11 @@ def get_unpacked_slice(self, index): post_index = index index //= 8 + return index, post_index + + def get_unpacked_slice(self, index): + # Numpy indexing is handled primarily in https://github.com/numpy/numpy/blob/maintenance/1.26.x/numpy/core/src/multiarray/mapping.c#L1435 + index, post_index = self.get_index_parameters(index) if post_index is NotImplemented: raise NotImplementedError(f"The following indexing scheme is not supported:\n{index}\n" "If you believe this scheme should be supported, " @@ -510,7 +526,24 @@ def __getitem__(self, item): return self.get_unpacked_slice(item) def set_unpacked_slice(self, index, value): - pass + assert not isinstance(value, GF2BP) + + packed_index, post_index = self.get_index_parameters(index) + + packed = self.view(np.ndarray)[packed_index] + if np.isscalar(packed): + packed = GF2BP([packed], self._axis_count).view(np.ndarray) + if packed.ndim == 1 and self.ndim > 1: + packed = packed[:, None] + + unpacked = np.unpackbits(packed, axis=-1, count=self._axis_count) + unpacked[post_index] = value + repacked = np.packbits(unpacked.view(np.ndarray), axis=-1) + + self.view(np.ndarray)[packed_index] = repacked[packed_index] + + def __setitem__(self, item, value): + self.set_unpacked_slice(item, value) GF2._default_ufunc_mode = "jit-calculate" diff --git a/tests/fields/test_bitpacked.py b/tests/fields/test_bitpacked.py index b11cc2336..d304bed11 100644 --- a/tests/fields/test_bitpacked.py +++ b/tests/fields/test_bitpacked.py @@ -55,11 +55,88 @@ def test_galois_array_indexing(): sub_matrix = arr_2d[np.ix_(row_indices, col_indices)] assert np.array_equal(sub_matrix, GF([[1, 0], [0, 1]])) + # TODO: Do we want to support this function? # 11. Indexing with np.take - taken = np.take(arr, [0, 2]) - assert np.array_equal(taken, GF([1, 1])) + # taken = np.take(arr, [0, 2]) + # assert np.array_equal(taken, GF([1, 1])) print("All tests passed.") + +def test_galois_array_setting(): + # Define a Galois field array + GF = galois.GF(2) + arr = GF([1, 0, 1, 1]) + arr = np.packbits(arr) + + # 1. Basic Indexing + arr[0] = GF(0) + assert arr[0] == GF(0) + + # 2. Negative Indexing + arr[-1] = GF(0) + assert arr[-1] == GF(0) + + # 3. Slicing + arr[1:3] = GF([1, 0]) + assert np.array_equal(arr, np.packbits(GF([0, 1, 0, 0]))) + + # 4. Multidimensional Indexing + arr_2d = GF([[1, 0], [0, 1]]) + arr_2d = np.packbits(arr_2d) + arr_2d[0, 1] = GF(1) + assert arr_2d[0, 1] == GF(1) + + arr_2d[:, 1] = GF([0, 0]) + assert np.array_equal(arr_2d[:, 1], GF([0, 0])) + + # 5. Boolean Indexing + arr = GF([1, 0, 1, 1]) + arr = np.packbits(arr) + mask = np.array([True, False, True, False]) + arr[mask] = GF(0) + assert np.array_equal(arr, np.packbits(GF([0, 0, 0, 1]))) + + # 6. Fancy Indexing + arr = GF([1, 0, 1, 1]) + arr = np.packbits(arr) + indices = [0, 2, 3] + arr[indices] = GF([0, 0, 0]) + assert np.array_equal(arr, np.packbits(GF([0, 0, 0, 0]))) + + # 7. Ellipsis + arr_3d = GF(np.random.randint(0, 2, (2, 3, 4))) + arr_3d = np.packbits(arr_3d) + arr_3d[0, ..., 1] = GF(1) + assert np.array_equal(arr_3d[0, :, 1], GF([1, 1, 1])) + + # 8. Indexing with slice objects + arr = GF([1, 0, 1, 1]) + arr = np.packbits(arr) + s = slice(1, 3) + arr[s] = GF([0, 0]) + assert np.array_equal(arr, np.packbits(GF([1, 0, 0, 1]))) + + # 9. Using np.newaxis (reshaped array assignment) + arr = GF([1, 0, 1, 1]) + arr = np.packbits(arr) + reshaped = arr[:, np.newaxis] # should this be using arr's data (as would be the case without packbits) or a new array? + reshaped = np.packbits(reshaped) + reshaped[:, 0] = GF([0, 0, 0, 0]) + # assert np.array_equal(arr, np.packbits(GF([0, 0, 0, 0]))) + assert np.array_equal(reshaped, np.packbits(GF([[0], [0], [0], [0]]))) + + # 10. Indexing with np.ix_ + arr_2d = GF([[1, 0], [0, 1]]) + arr_2d = np.packbits(arr_2d) + row_indices = np.array([0, 1]) + col_indices = np.array([0, 1]) + arr_2d[np.ix_(row_indices, col_indices)] = GF([[0, 0], [0, 0]]) + assert np.array_equal(arr_2d, np.packbits(GF([[0, 0], [0, 0]]))) + + print("All set-indexing tests passed.") + + if __name__ == "__main__": test_galois_array_indexing() + test_galois_array_setting() From 6c1d62fd4c3a6ca396ce5591de263cbb9ee4772c Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Tue, 24 Dec 2024 12:26:38 -0800 Subject: [PATCH 13/21] Get inv working --- src/galois/_fields/_gf2.py | 51 +++++++++++++++++++++++++++++++------- 1 file changed, 42 insertions(+), 9 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 0a9c8ec05..418a07975 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -226,6 +226,7 @@ class concatenate_bitpacked(Function): """Concatenates matrices together""" def __call__(self, arrays: list[Array], **kwargs): + # TODO: Should we only unpack arrays that have column counts %8 != 0 and concatenate those first? unpacked_arrays = [] for array in arrays: verify_isinstance(array, self.field) @@ -433,6 +434,20 @@ def __init__( ): pass + def __array_finalize__(self, obj): + """ + A NumPy dunder method that is called after "new", "view", or "new from template". It is used here to ensure + that view casting to a Galois field array has the appropriate dtype and that the values are in the field. + """ + super().__array_finalize__(obj) + + # Pass along _axis_count if it has been set already + try: + self._axis_count = obj._axis_count + except AttributeError: + pass + + @classmethod def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: r""" @@ -459,8 +474,20 @@ def get_index_parameters(self, index): post_index = (slice(None), col_index) index = (row_index, col_index // 8) elif isinstance(col_index, slice): - post_index = (slice(None), col_index) - col_index = slice(col_index.start // 8, max(col_index.step // 8, 1), max(col_index.stop // 8, 1)) + if isinstance(row_index, np.ndarray): + post_index = np.array(range(len(row_index))).reshape(row_index.shape), col_index + elif isinstance(row_index, int): + post_index = 0, col_index + elif isinstance(row_index, slice): + post_index = slice(0 if row_index.start is not None else row_index.start, + (row_index.stop - row_index.start) // row_index.step if row_index.stop is not None else row_index.stop, + 1 if row_index.step is not None else row_index.step), col_index + else: + post_index = list(range(len(row_index))), col_index + + col_index = slice(col_index.start // 8 if col_index.start is not None else col_index.start, + max(col_index.stop // 8, 1) if col_index.stop is not None else col_index.stop, + max(col_index.step // 8, 1) if col_index.step is not None else col_index.step) index = (row_index, col_index) elif isinstance(col_index, (Sequence, np.ndarray)): if isinstance(row_index, np.ndarray): @@ -483,8 +510,8 @@ def get_index_parameters(self, index): elif isinstance(index, tuple) and any(isinstance(x, slice) for x in index): post_index = index[1:] axis_adjustment = (slice(index.start // 8 if index.start is not None else index.start, - max(index.step // 8, 1) if index.step is not None else index.step, - max(index.stop // 8, 1) if index.stop is not None else index.stop) + max(index.stop // 8, 1) if index.stop is not None else index.stop, + max(index.step // 8, 1) if index.step is not None else index.step) if isinstance(index[-1], slice) else (index[-1] // 8,)) index = index[:-1] + axis_adjustment elif isinstance(index, slice): @@ -495,8 +522,8 @@ def get_index_parameters(self, index): # Array is 1-D, so we need to adjust post_index = index index = slice(index.start // 8 if index.start is not None else index.start, - max(index.step // 8, 1) if index.step is not None else index.step, - max(index.stop // 8, 1) if index.stop is not None else index.stop) + max(index.stop // 8, 1) if index.stop is not None else index.stop, + max(index.step // 8, 1) if index.step is not None else index.step) elif isinstance(index, int): post_index = index index //= 8 @@ -505,7 +532,7 @@ def get_index_parameters(self, index): def get_unpacked_slice(self, index): # Numpy indexing is handled primarily in https://github.com/numpy/numpy/blob/maintenance/1.26.x/numpy/core/src/multiarray/mapping.c#L1435 - index, post_index = self.get_index_parameters(index) + packed_index, post_index = self.get_index_parameters(index) if post_index is NotImplemented: raise NotImplementedError(f"The following indexing scheme is not supported:\n{index}\n" "If you believe this scheme should be supported, " @@ -514,7 +541,7 @@ def get_unpacked_slice(self, index): "`array = array.view(np.ndarray)` and then call the function." ) - packed = self.view(np.ndarray)[index] + packed = self.view(np.ndarray)[packed_index] if np.isscalar(packed): packed = GF2BP([packed], self._axis_count).view(np.ndarray) if packed.ndim == 1 and self.ndim > 1: @@ -539,8 +566,14 @@ def set_unpacked_slice(self, index, value): unpacked = np.unpackbits(packed, axis=-1, count=self._axis_count) unpacked[post_index] = value repacked = np.packbits(unpacked.view(np.ndarray), axis=-1) + if isinstance(packed_index[0], np.ndarray): + repacked_index = np.array(range(len(packed_index[0]))).reshape(packed_index.shape), packed_index[1] + elif isinstance(packed_index[0], int): + repacked_index = 0, packed_index[1] + else: + repacked_index = list(range(len(packed_index[0]))), packed_index[1] - self.view(np.ndarray)[packed_index] = repacked[packed_index] + self.view(np.ndarray)[packed_index] = repacked[repacked_index] def __setitem__(self, item, value): self.set_unpacked_slice(item, value) From f2b286a55992fa3f26a6f729c48d3adbf149128d Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Tue, 24 Dec 2024 13:07:39 -0800 Subject: [PATCH 14/21] Fix up repacked_index --- src/galois/_fields/_gf2.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 418a07975..9dac69f8a 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -566,12 +566,20 @@ def set_unpacked_slice(self, index, value): unpacked = np.unpackbits(packed, axis=-1, count=self._axis_count) unpacked[post_index] = value repacked = np.packbits(unpacked.view(np.ndarray), axis=-1) - if isinstance(packed_index[0], np.ndarray): - repacked_index = np.array(range(len(packed_index[0]))).reshape(packed_index.shape), packed_index[1] - elif isinstance(packed_index[0], int): - repacked_index = 0, packed_index[1] + if isinstance(packed_index, int) or isinstance(packed_index, slice) or len(packed_index) != 2: + repacked_index = packed_index else: - repacked_index = list(range(len(packed_index[0]))), packed_index[1] + packed_row_index, packed_col_index = packed_index + if isinstance(packed_row_index, np.ndarray): + repacked_index = np.array(range(len(packed_row_index))).reshape(packed_row_index.shape), packed_col_index + elif isinstance(packed_row_index, int): + repacked_index = 0, packed_col_index + elif isinstance(packed_row_index, slice): + repacked_index = slice(0 if packed_row_index.start is not None else packed_row_index.start, + (packed_row_index.stop - packed_row_index.start) // packed_row_index.step if packed_row_index.stop is not None else packed_row_index.stop, + 1 if packed_row_index.step is not None else packed_row_index.step), packed_col_index + else: + repacked_index = list(range(len(packed_row_index))), packed_col_index self.view(np.ndarray)[packed_index] = repacked[repacked_index] From 913b95df05500d9fad72ff184f72498df4c40ff0 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Thu, 26 Dec 2024 15:15:47 -0800 Subject: [PATCH 15/21] Implement bit-packed outer product --- src/galois/_fields/_gf2.py | 33 +++++++++++++++++++++++++++++++-- 1 file changed, 31 insertions(+), 2 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 9dac69f8a..638ba3be6 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -4,6 +4,8 @@ from __future__ import annotations +import operator +from functools import reduce from typing import Sequence import numpy as np @@ -163,8 +165,30 @@ class multiply_ufunc_bitpacked(multiply_ufunc): """ def __call__(self, ufunc, method, inputs, kwargs, meta): - output = super().__call__(ufunc, method, inputs, kwargs, meta) - output._axis_count = max(i._axis_count for i in inputs) + is_outer_product = method == "outer" + if is_outer_product and np.all(len(i.unpacked_shape) == 1 for i in inputs): + result_shape = reduce(operator.add, (i.unpacked_shape for i in inputs)) + else: + result_shape = np.broadcast_shapes(*(i.unpacked_shape for i in inputs)) + + if is_outer_product and len(inputs) == 2: + a = np.unpackbits(inputs[0]) + output = a[:, np.newaxis].view(np.ndarray) * inputs[1].view(np.ndarray) + else: + output = super().__call__(ufunc, method, inputs, kwargs, meta) + + assert len(output.shape) == len(result_shape) + # output = output.view(np.ndarray) + # if output.shape != result_shape: + # for axis, shape in enumerate(zip(output.shape, result_shape)): + # if axis == len(result_shape) - 1: + # # The last axis remains packed + # break + # + # if shape[0] != shape[1]: + # output = np.unpackbits(output, axis=axis, count=shape[1]) + output = self.field._view(output) + output._axis_count = result_shape[-1] return output @@ -530,6 +554,11 @@ def get_index_parameters(self, index): return index, post_index + # TODO: Should this be the default shape returned and a cast (i.e. .view(np.ndarray) is needed to get the packed shape? + @property + def unpacked_shape(self): + return self.shape[:-1] + (self._axis_count,) + def get_unpacked_slice(self, index): # Numpy indexing is handled primarily in https://github.com/numpy/numpy/blob/maintenance/1.26.x/numpy/core/src/multiarray/mapping.c#L1435 packed_index, post_index = self.get_index_parameters(index) From 28267d1d33ccec74c5b051b90e612412ba3e894d Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Fri, 27 Dec 2024 12:34:35 -0800 Subject: [PATCH 16/21] Try normalizing indexing --- src/galois/_fields/_gf2.py | 291 +++++++++++++++++++++++++------------ 1 file changed, 200 insertions(+), 91 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 638ba3be6..c524ab3d5 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -6,7 +6,7 @@ import operator from functools import reduce -from typing import Sequence +from typing import Sequence, Final import numpy as np from typing_extensions import Literal, Self, Optional @@ -275,6 +275,7 @@ def __call__(self, A: Array) -> Array: # Concatenate A and I to get the matrix AI = [A | I] AI = np.concatenate((A, I), axis=-1) + AI[0] # Perform Gaussian elimination to get the reduced row echelon form AI_rre = [I | A^-1] AI_rre, _ = row_reduce_jit(self.field)(AI, ncols=n) @@ -419,6 +420,7 @@ class factory :func:`~galois.GF`. Group: galois-fields """ + BIT_WIDTH = 8 def __new__( cls, @@ -489,70 +491,181 @@ def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: array = GF2.Identity(size, dtype=dtype) return np.packbits(array) + @staticmethod + def _normalize_indexing_to_tuple(index): + """ + Normalize indexing into a tuple of slices, integers, and ellipses. + + Args: + index: The indexing expression (int, slice, list, etc.). + + Returns: + A tuple of slices, integers, and ellipses. + """ + if isinstance(index, int): + return (index,) + elif isinstance(index, slice): + return (index,) + elif isinstance(index, list): + # Lists cannot be directly converted to slices, so leave as-is + return (index,) + elif isinstance(index, tuple): + normalized = [] + for i in index: + if isinstance(i, (int, slice, np.ndarray)) or i is Ellipsis or i is np.newaxis: + normalized.append(i) + else: + raise TypeError(f"Unsupported index type: {type(i)}") + return tuple(normalized) + elif index is Ellipsis: + return (Ellipsis,) + elif isinstance(index, (Sequence, np.ndarray)): + if index.dtype == bool: + # Boolean mask can remain as-is + return (index,) + else: + raise TypeError("Cannot normalize NumPy arrays directly into slices.") + else: + raise TypeError(f"Unsupported index type: {type(index)}") + def get_index_parameters(self, index): - post_index = NotImplemented - if isinstance(index, (Sequence, np.ndarray)): - if len(index) == 2: - row_index, col_index = index - if isinstance(col_index, int): - post_index = (slice(None), col_index) - index = (row_index, col_index // 8) - elif isinstance(col_index, slice): - if isinstance(row_index, np.ndarray): - post_index = np.array(range(len(row_index))).reshape(row_index.shape), col_index - elif isinstance(row_index, int): - post_index = 0, col_index - elif isinstance(row_index, slice): - post_index = slice(0 if row_index.start is not None else row_index.start, - (row_index.stop - row_index.start) // row_index.step if row_index.stop is not None else row_index.stop, - 1 if row_index.step is not None else row_index.step), col_index + # Convert the many different forms of indexing into a normalized form of slices + # normalized_index = tuple() + # if not isinstance(index, tuple): + # index = (index,) + + normalized_index = self._normalize_indexing_to_tuple(index) + + assert isinstance(normalized_index, tuple) + + bit_width: Final[int] = self.BIT_WIDTH + packed_index = tuple() + packed_stop = -1 + unpacked_index = tuple() + unpacked_start = -1 + for axis, i in enumerate(normalized_index): + is_unpacked_axis = axis != self.ndim - 1 # only the last axis is packed + if isinstance(i, int): + if is_unpacked_axis and self.ndim > 1: + packed_index += (i,) + unpacked_index = (0,) + else: + packed_index += (i // bit_width,) + unpacked_index += (i,) + elif isinstance(i, slice): + if is_unpacked_axis: + packed_index += (i,) + packed_stop = -1 + unpacked_index += (slice(0 if i.start is not None else i.start, + (i.stop - i.start) // i.step if i.stop is not None else i.stop, + 1 if i.step is not None else i.step),) + unpacked_start = -2 + else: + packed_index += (slice(i.start // bit_width if i.start is not None else i.start, + max(i.stop // bit_width, 1) if i.stop is not None else i.stop, + max(i.step // bit_width, 1) if i.step is not None else i.step),) + unpacked_index += (i,) + elif isinstance(i, (Sequence, np.ndarray)): + if is_unpacked_axis: + packed_index += (i,) + if isinstance(i, np.ndarray): + unpacked_index += (np.array(range(len(i))).reshape(i.shape),) + if i.ndim > 1: + packed_stop = None + unpacked_start = None else: - post_index = list(range(len(row_index))), col_index - - col_index = slice(col_index.start // 8 if col_index.start is not None else col_index.start, - max(col_index.stop // 8, 1) if col_index.stop is not None else col_index.stop, - max(col_index.step // 8, 1) if col_index.step is not None else col_index.step) - index = (row_index, col_index) - elif isinstance(col_index, (Sequence, np.ndarray)): - if isinstance(row_index, np.ndarray): - post_index = np.array(range(len(row_index))).reshape(row_index.shape), col_index + unpacked_index += np.array(range(len(i))) + else: + if all(isinstance(x, int) for x in index): + packed_index += ((s // bit_width for s in i),) else: - post_index = list(range(len(row_index))), col_index - col_index = tuple(s // 8 for s in col_index) - index = (row_index, col_index) - elif col_index is None: # new axis - post_index = (slice(None), None) - index = (row_index,) - elif ((isinstance(index, np.ndarray) and index.ndim == 1) or - (isinstance(index, list) and all(isinstance(x, int) for x in index))): - post_index = index - index = list(range((len(index) // 8) + 1)) - elif isinstance(index, tuple) and any(x is Ellipsis for x in index): - post_index = index[1:] - axis_adjustment = (slice(None),) if index[-1] is Ellipsis else (index[-1] // 8,) - index = index[:-1] + axis_adjustment - elif isinstance(index, tuple) and any(isinstance(x, slice) for x in index): - post_index = index[1:] - axis_adjustment = (slice(index.start // 8 if index.start is not None else index.start, - max(index.stop // 8, 1) if index.stop is not None else index.stop, - max(index.step // 8, 1) if index.step is not None else index.step) - if isinstance(index[-1], slice) else (index[-1] // 8,)) - index = index[:-1] + axis_adjustment - elif isinstance(index, slice): - if self.ndim > 1: - # Rows aren't packed, so we can index normally - post_index = slice(None) - if len(self.shape) == 1: - # Array is 1-D, so we need to adjust - post_index = index - index = slice(index.start // 8 if index.start is not None else index.start, - max(index.stop // 8, 1) if index.stop is not None else index.stop, - max(index.step // 8, 1) if index.step is not None else index.step) - elif isinstance(index, int): - post_index = index - index //= 8 - - return index, post_index + packed_index += (i // bit_width,) + + if isinstance(i, np.ndarray) and i.ndim > 1: + packed_stop = None + unpacked_start = None + + unpacked_index += (i,) + elif i is np.newaxis: + packed_index += (i,) + packed_stop = -1 + unpacked_index += (i,) + unpacked_start = -2 + elif i is Ellipsis: + packed_index += (i,) + packed_stop = -1 + unpacked_index += (i,) + unpacked_start = -2 + else: + raise NotImplementedError(f"The following indexing scheme is not supported:\n{index}\n" + "If you believe this scheme should be supported, " + "please submit a GitHub issue at https://github.com/mhostetter/galois/issues.\n\n" + "If you'd like to perform this operation on the data, you should first call " + "`array = array.view(np.ndarray)` and then call the function." + ) + + # if len(index) == 2: + # row_index, col_index = index + # if isinstance(col_index, int): + # post_index = (slice(None), col_index) + # index = (row_index, col_index // 8) + # elif isinstance(col_index, slice): + # if isinstance(row_index, np.ndarray): + # post_index = np.array(range(len(row_index))).reshape(row_index.shape), col_index + # elif isinstance(row_index, int): + # post_index = 0, col_index + # elif isinstance(row_index, slice): + # post_index = slice(0 if row_index.start is not None else row_index.start, + # (row_index.stop - row_index.start) // row_index.step if row_index.stop is not None else row_index.stop, + # 1 if row_index.step is not None else row_index.step), col_index + # else: + # post_index = list(range(len(row_index))), col_index + # + # col_index = slice(col_index.start // 8 if col_index.start is not None else col_index.start, + # max(col_index.stop // 8, 1) if col_index.stop is not None else col_index.stop, + # max(col_index.step // 8, 1) if col_index.step is not None else col_index.step) + # index = (row_index, col_index) + # elif isinstance(col_index, (Sequence, np.ndarray)): + # if isinstance(row_index, np.ndarray): + # post_index = np.array(range(len(row_index))).reshape(row_index.shape), col_index + # else: + # post_index = list(range(len(row_index))), col_index + # col_index = tuple(s // 8 for s in col_index) + # index = (row_index, col_index) + # elif col_index is None: # new axis + # post_index = (slice(None), None) + # index = (row_index,) + # elif ((isinstance(index, np.ndarray) and index.ndim == 1) or + # (isinstance(index, list) and all(isinstance(x, int) for x in index))): + # post_index = index + # index = list(range((len(index) // 8) + 1)) + # elif isinstance(index, tuple) and any(x is Ellipsis for x in index): + # post_index = index[1:] + # axis_adjustment = (slice(None),) if index[-1] is Ellipsis else (index[-1] // 8,) + # index = index[:-1] + axis_adjustment + # elif isinstance(index, tuple) and any(isinstance(x, slice) for x in index): + # post_index = index[1:] + # axis_adjustment = (slice(index.start // 8 if index.start is not None else index.start, + # max(index.stop // 8, 1) if index.stop is not None else index.stop, + # max(index.step // 8, 1) if index.step is not None else index.step) + # if isinstance(index[-1], slice) else (index[-1] // 8,)) + # index = index[:-1] + axis_adjustment + # elif isinstance(index, slice): + # if self.ndim > 1: + # # Rows aren't packed, so we can index normally + # post_index = slice(None) + # if len(self.shape) == 1: + # # Array is 1-D, so we need to adjust + # post_index = index + # index = slice(index.start // 8 if index.start is not None else index.start, + # max(index.stop // 8, 1) if index.stop is not None else index.stop, + # max(index.step // 8, 1) if index.step is not None else index.step) + # elif isinstance(index, int): + # post_index = index + # index //= 8 + + return (packed_index[:packed_stop] if packed_stop is not None else packed_index, + unpacked_index[unpacked_start:] if unpacked_start is not None else tuple()) # TODO: Should this be the default shape returned and a cast (i.e. .view(np.ndarray) is needed to get the packed shape? @property @@ -561,22 +674,17 @@ def unpacked_shape(self): def get_unpacked_slice(self, index): # Numpy indexing is handled primarily in https://github.com/numpy/numpy/blob/maintenance/1.26.x/numpy/core/src/multiarray/mapping.c#L1435 - packed_index, post_index = self.get_index_parameters(index) - if post_index is NotImplemented: - raise NotImplementedError(f"The following indexing scheme is not supported:\n{index}\n" - "If you believe this scheme should be supported, " - "please submit a GitHub issue at https://github.com/mhostetter/galois/issues.\n\n" - "If you'd like to perform this operation on the data, you should first call " - "`array = array.view(np.ndarray)` and then call the function." - ) + packed_index, unpacked_index = self.get_index_parameters(index) packed = self.view(np.ndarray)[packed_index] if np.isscalar(packed): packed = GF2BP([packed], self._axis_count).view(np.ndarray) - if packed.ndim == 1 and self.ndim > 1: - packed = packed[:, None] + # if packed.ndim == 1 and self.ndim > 1: # for some reason this is causing an issue when the column count > 8 + # packed = packed[:, None] unpacked = np.unpackbits(packed, axis=-1, count=self._axis_count) - return GF2._view(unpacked[post_index]) + # if unpacked.ndim == 1 and self.ndim > 1: # for some reason this is causing an issue when the column count > 8 + # unpacked = unpacked[None, :] + return GF2._view(unpacked[unpacked_index]) def __getitem__(self, item): return self.get_unpacked_slice(item) @@ -589,28 +697,29 @@ def set_unpacked_slice(self, index, value): packed = self.view(np.ndarray)[packed_index] if np.isscalar(packed): packed = GF2BP([packed], self._axis_count).view(np.ndarray) - if packed.ndim == 1 and self.ndim > 1: - packed = packed[:, None] + # if packed.ndim == 1 and self.ndim > 1: + # packed = packed[:, None] unpacked = np.unpackbits(packed, axis=-1, count=self._axis_count) unpacked[post_index] = value repacked = np.packbits(unpacked.view(np.ndarray), axis=-1) - if isinstance(packed_index, int) or isinstance(packed_index, slice) or len(packed_index) != 2: - repacked_index = packed_index - else: - packed_row_index, packed_col_index = packed_index - if isinstance(packed_row_index, np.ndarray): - repacked_index = np.array(range(len(packed_row_index))).reshape(packed_row_index.shape), packed_col_index - elif isinstance(packed_row_index, int): - repacked_index = 0, packed_col_index - elif isinstance(packed_row_index, slice): - repacked_index = slice(0 if packed_row_index.start is not None else packed_row_index.start, - (packed_row_index.stop - packed_row_index.start) // packed_row_index.step if packed_row_index.stop is not None else packed_row_index.stop, - 1 if packed_row_index.step is not None else packed_row_index.step), packed_col_index - else: - repacked_index = list(range(len(packed_row_index))), packed_col_index - - self.view(np.ndarray)[packed_index] = repacked[repacked_index] + repacked_index = packed_index + # if isinstance(packed_index, int) or isinstance(packed_index, slice) or len(packed_index) != 2: + # repacked_index = packed_index + # else: + # packed_row_index, packed_col_index = packed_index + # if isinstance(packed_row_index, np.ndarray): + # repacked_index = np.array(range(len(packed_row_index))).reshape(packed_row_index.shape), packed_col_index + # elif isinstance(packed_row_index, int): + # repacked_index = 0, packed_col_index + # elif isinstance(packed_row_index, slice): + # repacked_index = slice(0 if packed_row_index.start is not None else packed_row_index.start, + # (packed_row_index.stop - packed_row_index.start) // packed_row_index.step if packed_row_index.stop is not None else packed_row_index.stop, + # 1 if packed_row_index.step is not None else packed_row_index.step), packed_col_index + # else: + # repacked_index = list(range(len(packed_row_index))), packed_col_index + + self.view(np.ndarray)[packed_index] = repacked def __setitem__(self, item, value): self.set_unpacked_slice(item, value) From 7d44824662cdbfba9ba201ddc7341e5ddb251839 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Tue, 31 Dec 2024 09:35:38 -0800 Subject: [PATCH 17/21] Clean up indexing rules --- src/galois/_fields/_gf2.py | 223 +++++++++++++------------------------ 1 file changed, 79 insertions(+), 144 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index c524ab3d5..2d76e17d6 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -492,74 +492,78 @@ def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: return np.packbits(array) @staticmethod - def _normalize_indexing_to_tuple(index): + def _normalize_indexing_to_tuple(index, ndim): """ - Normalize indexing into a tuple of slices, integers, and ellipses. + Normalize indexing into a tuple of slices, integers, and/or new axes. + NOTE: Ellipsis indexing is converted to slice indexing. Args: index: The indexing expression (int, slice, list, etc.). + ndim: The number of dimensions of the array being indexed into. Returns: - A tuple of slices, integers, and ellipses. + A tuple of integers, slices, and/or new axes. """ if isinstance(index, int): return (index,) elif isinstance(index, slice): return (index,) elif isinstance(index, list): - # Lists cannot be directly converted to slices, so leave as-is + # Lists cannot be directly converted to slices, so leave as-is. + return (index,) + elif index is np.newaxis: return (index,) elif isinstance(index, tuple): normalized = [] - for i in index: - if isinstance(i, (int, slice, np.ndarray)) or i is Ellipsis or i is np.newaxis: - normalized.append(i) - else: - raise TypeError(f"Unsupported index type: {type(i)}") + + if any(i is Ellipsis for i in index): + num_explicit_dims = sum(1 for i in index if i is not Ellipsis) + for i in index: + if i is Ellipsis: + normalized.extend([slice(None)] * (ndim - num_explicit_dims)) + else: + normalized.append(i) + else: + for i in index: + normalized.extend(GF2BP._normalize_indexing_to_tuple(i, ndim)) + return tuple(normalized) - elif index is Ellipsis: - return (Ellipsis,) elif isinstance(index, (Sequence, np.ndarray)): - if index.dtype == bool: - # Boolean mask can remain as-is - return (index,) - else: - raise TypeError("Cannot normalize NumPy arrays directly into slices.") + # Boolean mask or fancy indexing can remain as-is + return (index,) else: - raise TypeError(f"Unsupported index type: {type(index)}") + raise TypeError(f"Unsupported indexing type: {type(index)}") def get_index_parameters(self, index): - # Convert the many different forms of indexing into a normalized form of slices - # normalized_index = tuple() - # if not isinstance(index, tuple): - # index = (index,) - - normalized_index = self._normalize_indexing_to_tuple(index) + normalized_index = self._normalize_indexing_to_tuple(index, self.ndim) assert isinstance(normalized_index, tuple) bit_width: Final[int] = self.BIT_WIDTH packed_index = tuple() - packed_stop = -1 unpacked_index = tuple() - unpacked_start = -1 + shape = tuple() + axes_in_index = len(normalized_index) for axis, i in enumerate(normalized_index): is_unpacked_axis = axis != self.ndim - 1 # only the last axis is packed if isinstance(i, int): + shape += (1,) if is_unpacked_axis and self.ndim > 1: packed_index += (i,) - unpacked_index = (0,) + + # If we have multidimensional indexing, then we will need to re-index the same way after reshaping + if axes_in_index > 1: + unpacked_index = (0,) else: packed_index += (i // bit_width,) unpacked_index += (i,) elif isinstance(i, slice): + shape += (self.shape[axis],) if is_unpacked_axis: packed_index += (i,) - packed_stop = -1 unpacked_index += (slice(0 if i.start is not None else i.start, (i.stop - i.start) // i.step if i.stop is not None else i.stop, 1 if i.step is not None else i.step),) - unpacked_start = -2 else: packed_index += (slice(i.start // bit_width if i.start is not None else i.start, max(i.stop // bit_width, 1) if i.stop is not None else i.stop, @@ -568,34 +572,36 @@ def get_index_parameters(self, index): elif isinstance(i, (Sequence, np.ndarray)): if is_unpacked_axis: packed_index += (i,) - if isinstance(i, np.ndarray): - unpacked_index += (np.array(range(len(i))).reshape(i.shape),) - if i.ndim > 1: - packed_stop = None - unpacked_start = None - else: - unpacked_index += np.array(range(len(i))) else: - if all(isinstance(x, int) for x in index): - packed_index += ((s // bit_width for s in i),) + if isinstance(index, np.ndarray) and index.dtype == np.bool: + mask_packed = [False] * self.shape[axis] + for j, value in enumerate(i): + mask_packed[j // bit_width] |= True + packed_index = mask_packed + unpacked_index = i + shape += (max(sum(i) // bit_width, 1),) else: - packed_index += (i // bit_width,) - - if isinstance(i, np.ndarray) and i.ndim > 1: - packed_stop = None - unpacked_start = None - - unpacked_index += (i,) + # adjust indexing for this packed axis + data = np.array([s // bit_width for s in i], dtype=self.dtype) + # remove duplicate entries, including for nested arrays + if data.ndim > 1: + rows = [] + for j, row_data in enumerate(data): + _, unique_indices = np.unique(row_data, return_index=True) + # maintain the original order + rows.append(row_data[np.sort(unique_indices)]) + data = np.vstack(rows) + else: + _, unique_indices = np.unique(data, return_index=True) + # maintain the original order + data = data[np.sort(unique_indices)] + + packed_index += (data,) + if axes_in_index == 1: + unpacked_index += (i,) elif i is np.newaxis: packed_index += (i,) - packed_stop = -1 unpacked_index += (i,) - unpacked_start = -2 - elif i is Ellipsis: - packed_index += (i,) - packed_stop = -1 - unpacked_index += (i,) - unpacked_start = -2 else: raise NotImplementedError(f"The following indexing scheme is not supported:\n{index}\n" "If you believe this scheme should be supported, " @@ -604,131 +610,60 @@ def get_index_parameters(self, index): "`array = array.view(np.ndarray)` and then call the function." ) - # if len(index) == 2: - # row_index, col_index = index - # if isinstance(col_index, int): - # post_index = (slice(None), col_index) - # index = (row_index, col_index // 8) - # elif isinstance(col_index, slice): - # if isinstance(row_index, np.ndarray): - # post_index = np.array(range(len(row_index))).reshape(row_index.shape), col_index - # elif isinstance(row_index, int): - # post_index = 0, col_index - # elif isinstance(row_index, slice): - # post_index = slice(0 if row_index.start is not None else row_index.start, - # (row_index.stop - row_index.start) // row_index.step if row_index.stop is not None else row_index.stop, - # 1 if row_index.step is not None else row_index.step), col_index - # else: - # post_index = list(range(len(row_index))), col_index - # - # col_index = slice(col_index.start // 8 if col_index.start is not None else col_index.start, - # max(col_index.stop // 8, 1) if col_index.stop is not None else col_index.stop, - # max(col_index.step // 8, 1) if col_index.step is not None else col_index.step) - # index = (row_index, col_index) - # elif isinstance(col_index, (Sequence, np.ndarray)): - # if isinstance(row_index, np.ndarray): - # post_index = np.array(range(len(row_index))).reshape(row_index.shape), col_index - # else: - # post_index = list(range(len(row_index))), col_index - # col_index = tuple(s // 8 for s in col_index) - # index = (row_index, col_index) - # elif col_index is None: # new axis - # post_index = (slice(None), None) - # index = (row_index,) - # elif ((isinstance(index, np.ndarray) and index.ndim == 1) or - # (isinstance(index, list) and all(isinstance(x, int) for x in index))): - # post_index = index - # index = list(range((len(index) // 8) + 1)) - # elif isinstance(index, tuple) and any(x is Ellipsis for x in index): - # post_index = index[1:] - # axis_adjustment = (slice(None),) if index[-1] is Ellipsis else (index[-1] // 8,) - # index = index[:-1] + axis_adjustment - # elif isinstance(index, tuple) and any(isinstance(x, slice) for x in index): - # post_index = index[1:] - # axis_adjustment = (slice(index.start // 8 if index.start is not None else index.start, - # max(index.stop // 8, 1) if index.stop is not None else index.stop, - # max(index.step // 8, 1) if index.step is not None else index.step) - # if isinstance(index[-1], slice) else (index[-1] // 8,)) - # index = index[:-1] + axis_adjustment - # elif isinstance(index, slice): - # if self.ndim > 1: - # # Rows aren't packed, so we can index normally - # post_index = slice(None) - # if len(self.shape) == 1: - # # Array is 1-D, so we need to adjust - # post_index = index - # index = slice(index.start // 8 if index.start is not None else index.start, - # max(index.stop // 8, 1) if index.stop is not None else index.stop, - # max(index.step // 8, 1) if index.step is not None else index.step) - # elif isinstance(index, int): - # post_index = index - # index //= 8 - - return (packed_index[:packed_stop] if packed_stop is not None else packed_index, - unpacked_index[unpacked_start:] if unpacked_start is not None else tuple()) + return packed_index, unpacked_index, shape # TODO: Should this be the default shape returned and a cast (i.e. .view(np.ndarray) is needed to get the packed shape? @property def unpacked_shape(self): return self.shape[:-1] + (self._axis_count,) - def get_unpacked_slice(self, index): + def get_unpacked_value(self, index): # Numpy indexing is handled primarily in https://github.com/numpy/numpy/blob/maintenance/1.26.x/numpy/core/src/multiarray/mapping.c#L1435 - packed_index, unpacked_index = self.get_index_parameters(index) + packed_index, unpacked_index, shape = self.get_index_parameters(index) packed = self.view(np.ndarray)[packed_index] + if np.isscalar(packed): packed = GF2BP([packed], self._axis_count).view(np.ndarray) - # if packed.ndim == 1 and self.ndim > 1: # for some reason this is causing an issue when the column count > 8 - # packed = packed[:, None] + + if len(shape) > 0: + packed = packed.reshape(shape) + unpacked = np.unpackbits(packed, axis=-1, count=self._axis_count) - # if unpacked.ndim == 1 and self.ndim > 1: # for some reason this is causing an issue when the column count > 8 - # unpacked = unpacked[None, :] return GF2._view(unpacked[unpacked_index]) def __getitem__(self, item): - return self.get_unpacked_slice(item) + return self.get_unpacked_value(item) - def set_unpacked_slice(self, index, value): + def set_unpacked_value(self, index, value): assert not isinstance(value, GF2BP) - packed_index, post_index = self.get_index_parameters(index) + packed_index, unpacked_index, shape = self.get_index_parameters(index) packed = self.view(np.ndarray)[packed_index] + original_packed_shape = packed.shape + if np.isscalar(packed): packed = GF2BP([packed], self._axis_count).view(np.ndarray) - # if packed.ndim == 1 and self.ndim > 1: - # packed = packed[:, None] + + if len(shape) > 0: + packed = packed.reshape(shape) unpacked = np.unpackbits(packed, axis=-1, count=self._axis_count) - unpacked[post_index] = value + unpacked[unpacked_index] = value repacked = np.packbits(unpacked.view(np.ndarray), axis=-1) - repacked_index = packed_index - # if isinstance(packed_index, int) or isinstance(packed_index, slice) or len(packed_index) != 2: - # repacked_index = packed_index - # else: - # packed_row_index, packed_col_index = packed_index - # if isinstance(packed_row_index, np.ndarray): - # repacked_index = np.array(range(len(packed_row_index))).reshape(packed_row_index.shape), packed_col_index - # elif isinstance(packed_row_index, int): - # repacked_index = 0, packed_col_index - # elif isinstance(packed_row_index, slice): - # repacked_index = slice(0 if packed_row_index.start is not None else packed_row_index.start, - # (packed_row_index.stop - packed_row_index.start) // packed_row_index.step if packed_row_index.stop is not None else packed_row_index.stop, - # 1 if packed_row_index.step is not None else packed_row_index.step), packed_col_index - # else: - # repacked_index = list(range(len(packed_row_index))), packed_col_index + repacked = repacked.reshape(original_packed_shape) self.view(np.ndarray)[packed_index] = repacked def __setitem__(self, item, value): - self.set_unpacked_slice(item, value) + self.set_unpacked_value(item, value) GF2._default_ufunc_mode = "jit-calculate" GF2._ufunc_modes = ["jit-calculate", "python-calculate"] GF2.compile("auto") -# GF2BP._default_ufunc_mode = "jit-calculate" -# GF2BP._ufunc_modes = ["jit-calculate", "python-calculate"] -# GF2BP.compile("auto") +GF2BP._default_ufunc_mode = "jit-calculate" +GF2BP._ufunc_modes = ["jit-calculate", "python-calculate"] +GF2BP.compile("auto") From 57bf589e8fc2536bda4bdb018e3e4273967d73b2 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Tue, 31 Dec 2024 11:19:59 -0800 Subject: [PATCH 18/21] Fix shape / slice bugs --- src/galois/_fields/_gf2.py | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 2d76e17d6..6aec3f163 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -547,31 +547,42 @@ def get_index_parameters(self, index): for axis, i in enumerate(normalized_index): is_unpacked_axis = axis != self.ndim - 1 # only the last axis is packed if isinstance(i, int): - shape += (1,) if is_unpacked_axis and self.ndim > 1: packed_index += (i,) # If we have multidimensional indexing, then we will need to re-index the same way after reshaping if axes_in_index > 1: unpacked_index = (0,) + shape += (1,) else: packed_index += (i // bit_width,) unpacked_index += (i,) + if axes_in_index > 1: + shape += (1,) + else: + shape += (self.shape[axis],) elif isinstance(i, slice): - shape += (self.shape[axis],) if is_unpacked_axis: packed_index += (i,) - unpacked_index += (slice(0 if i.start is not None else i.start, - (i.stop - i.start) // i.step if i.stop is not None else i.stop, - 1 if i.step is not None else i.step),) + # the packed index will already filter, so we just select everything after + unpacked_index += (slice(None),) else: packed_index += (slice(i.start // bit_width if i.start is not None else i.start, max(i.stop // bit_width, 1) if i.stop is not None else i.stop, max(i.step // bit_width, 1) if i.step is not None else i.step),) unpacked_index += (i,) + + packed_slice = packed_index[-1] + packed_slice = slice(0 if packed_slice.start is None else packed_slice.start, + self.shape[axis] if packed_slice.stop is None else packed_slice.stop, + 1 if packed_slice.step is None else packed_slice.step) + slice_size = max(0, (packed_slice.stop - packed_slice.start + packed_slice.step - 1) // packed_slice.step) + shape += (slice_size,) elif isinstance(i, (Sequence, np.ndarray)): if is_unpacked_axis: + shape += (len(i),) packed_index += (i,) + unpacked_index += (slice(None),) else: if isinstance(index, np.ndarray) and index.dtype == np.bool: mask_packed = [False] * self.shape[axis] @@ -597,6 +608,7 @@ def get_index_parameters(self, index): data = data[np.sort(unique_indices)] packed_index += (data,) + shape += (self.shape[axis],) if axes_in_index == 1: unpacked_index += (i,) elif i is np.newaxis: @@ -630,7 +642,11 @@ def get_unpacked_value(self, index): packed = packed.reshape(shape) unpacked = np.unpackbits(packed, axis=-1, count=self._axis_count) - return GF2._view(unpacked[unpacked_index]) + value = unpacked[unpacked_index] + if np.isscalar(value): + return GF2(value, dtype=self.dtype) + + return GF2._view(value) def __getitem__(self, item): return self.get_unpacked_value(item) From 8d2ad2659c0633bf4d867e0facbe0cb32c6f028c Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Thu, 9 Jan 2025 20:03:02 -0800 Subject: [PATCH 19/21] Normalize to positive indexing; .shape is now the unpacked shape; fix indexing bugs --- src/galois/_fields/_gf2.py | 126 ++++++++++++++++++++++++------------- 1 file changed, 81 insertions(+), 45 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 6aec3f163..04fc16d3e 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -6,6 +6,7 @@ import operator from functools import reduce +from math import ceil, floor from typing import Sequence, Final import numpy as np @@ -166,10 +167,10 @@ class multiply_ufunc_bitpacked(multiply_ufunc): def __call__(self, ufunc, method, inputs, kwargs, meta): is_outer_product = method == "outer" - if is_outer_product and np.all(len(i.unpacked_shape) == 1 for i in inputs): - result_shape = reduce(operator.add, (i.unpacked_shape for i in inputs)) + if is_outer_product and np.all(len(i.shape) == 1 for i in inputs): + result_shape = reduce(operator.add, (i.shape for i in inputs)) else: - result_shape = np.broadcast_shapes(*(i.unpacked_shape for i in inputs)) + result_shape = np.broadcast_shapes(*(i.shape for i in inputs)) if is_outer_product and len(inputs) == 2: a = np.unpackbits(inputs[0]) @@ -177,7 +178,7 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): else: output = super().__call__(ufunc, method, inputs, kwargs, meta) - assert len(output.shape) == len(result_shape) + assert len(output.view(np.ndarray).shape) == len(result_shape) # output = output.view(np.ndarray) # if output.shape != result_shape: # for axis, shape in enumerate(zip(output.shape, result_shape)): @@ -225,11 +226,13 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): b._axis_count = row_axis_count # Make sure the inner dimensions match (e.g. (M, N) x (N, P) -> (M, P)) - assert a.shape[-1] == b.shape[0] - if len(b.shape) == 1: - final_shape = (a.shape[0],) + a_packed_shape = a.view(np.ndarray).shape + b_packed_shape = b.view(np.ndarray).shape + assert a_packed_shape[-1] == b_packed_shape[0] + if len(b_packed_shape) == 1: + final_shape = (a_packed_shape[0],) else: - final_shape = (a.shape[0], b.shape[-1]) + final_shape = (a_packed_shape[0], b_packed_shape[-1]) if len(b.shape) == 1: # matrix-vector multiplication @@ -237,10 +240,11 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): else: # matrix-matrix multiplication output = GF2.Zeros(final_shape) - for i in range(b.shape[-1]): + for i in range(b_packed_shape[-1]): # TODO: Include alternate path for numpy < v2 # output[:, i] = np.bitwise_xor.reduce(np.unpackbits((a & b[:, i]).view(np.ndarray), axis=-1), axis=-1) - output[:, i] = np.bitwise_xor.reduce(np.bitwise_count((a & b[:, i]).view(np.ndarray)), axis=-1) % 2 + output[:, i] = np.bitwise_xor.reduce( + np.bitwise_count((a & b.view(np.ndarray)[:, i]).view(np.ndarray)), axis=-1) % 2 output = field._view(np.packbits(output.view(np.ndarray), axis=-1)) output._axis_count = final_shape[-1] @@ -275,7 +279,6 @@ def __call__(self, A: Array) -> Array: # Concatenate A and I to get the matrix AI = [A | I] AI = np.concatenate((A, I), axis=-1) - AI[0] # Perform Gaussian elimination to get the reduced row echelon form AI_rre = [I | A^-1] AI_rre, _ = row_reduce_jit(self.field)(AI, ncols=n) @@ -492,40 +495,67 @@ def Identity(cls, size: int, dtype: DTypeLike | None = None) -> Self: return np.packbits(array) @staticmethod - def _normalize_indexing_to_tuple(index, ndim): + def _normalize_indexing_to_tuple(index, shape, axis = 0): """ - Normalize indexing into a tuple of slices, integers, and/or new axes. + Normalize indexing into a tuple of positive-only slices, integers, and/or new axes. NOTE: Ellipsis indexing is converted to slice indexing. Args: index: The indexing expression (int, slice, list, etc.). - ndim: The number of dimensions of the array being indexed into. + shape: Tuple of integers representing the shape of the object being indexed. Returns: - A tuple of integers, slices, and/or new axes. + A tuple of positive integers, slices, and/or new axes. """ + if not isinstance(shape, tuple): + raise TypeError("Shape must be a tuple of integers.") + + ndim = len(shape) + if isinstance(index, int): + if index < 0: + index += shape[axis] return (index,) elif isinstance(index, slice): - return (index,) + start, stop, step = index.start, index.stop, index.step + step = step if step is not None else 1 + + if step > 0: + start = start if start is not None else 0 + stop = stop if stop is not None else shape[axis] + + # Adjust negative start/stop values + if start < 0: + start += shape[axis] + if stop < 0: + stop += shape[axis] + else: + start = start if start is not None else shape[axis] - 1 + stop = stop if stop is not None else -shape[axis] - 1 + + return (slice(start, stop, step),) elif isinstance(index, list): - # Lists cannot be directly converted to slices, so leave as-is. + for axis, i in enumerate(index): + if i < 0: + index[axis] += shape[axis] return (index,) elif index is np.newaxis: return (index,) elif isinstance(index, tuple): normalized = [] - if any(i is Ellipsis for i in index): - num_explicit_dims = sum(1 for i in index if i is not Ellipsis) - for i in index: - if i is Ellipsis: - normalized.extend([slice(None)] * (ndim - num_explicit_dims)) - else: - normalized.append(i) - else: - for i in index: - normalized.extend(GF2BP._normalize_indexing_to_tuple(i, ndim)) + num_explicit_dims = sum(1 for i in index if i is not Ellipsis) + for i in index: + if i is Ellipsis: + span = ndim - num_explicit_dims + expanded_dims = [slice(None)] * span + for e_axis, e in enumerate(expanded_dims): + expanded_dims[e_axis] = GF2BP._normalize_indexing_to_tuple(e, shape, axis)[0] + axis += 1 + normalized.extend(expanded_dims) + else: + normalized.extend(GF2BP._normalize_indexing_to_tuple(i, shape, axis)) + axis += 1 return tuple(normalized) elif isinstance(index, (Sequence, np.ndarray)): @@ -535,11 +565,12 @@ def _normalize_indexing_to_tuple(index, ndim): raise TypeError(f"Unsupported indexing type: {type(index)}") def get_index_parameters(self, index): - normalized_index = self._normalize_indexing_to_tuple(index, self.ndim) + normalized_index = self._normalize_indexing_to_tuple(index, self.shape) assert isinstance(normalized_index, tuple) bit_width: Final[int] = self.BIT_WIDTH + packed_shape = self.view(np.ndarray).shape packed_index = tuple() unpacked_index = tuple() shape = tuple() @@ -556,27 +587,32 @@ def get_index_parameters(self, index): shape += (1,) else: packed_index += (i // bit_width,) - unpacked_index += (i,) + unpacked_index += (i % bit_width,) if axes_in_index > 1: shape += (1,) else: - shape += (self.shape[axis],) + shape += (packed_shape[axis],) elif isinstance(i, slice): if is_unpacked_axis: packed_index += (i,) # the packed index will already filter, so we just select everything after unpacked_index += (slice(None),) else: - packed_index += (slice(i.start // bit_width if i.start is not None else i.start, - max(i.stop // bit_width, 1) if i.stop is not None else i.stop, - max(i.step // bit_width, 1) if i.step is not None else i.step),) - unpacked_index += (i,) + if i.step > 0: + packed_index += (slice(i.start // bit_width, + max(int(ceil(i.stop / bit_width)), 1), + max(i.step // bit_width, 1)),) + unpacked_index += (slice(i.start % bit_width, i.start % bit_width + i.stop - i.start, i.step),) + else: + packed_index += (slice(i.start // bit_width, + max(int(floor(i.stop / bit_width)), -packed_shape[axis] -1), + min(i.step // bit_width, -1)),) + unpacked_index += (slice(i.start % bit_width, i.start % bit_width + i.stop - i.start, i.step),) + packed_slice = packed_index[-1] - packed_slice = slice(0 if packed_slice.start is None else packed_slice.start, - self.shape[axis] if packed_slice.stop is None else packed_slice.stop, - 1 if packed_slice.step is None else packed_slice.step) - slice_size = max(0, (packed_slice.stop - packed_slice.start + packed_slice.step - 1) // packed_slice.step) + abs_step = abs(packed_slice.step) + slice_size = max(0, (packed_slice.stop - packed_slice.start + abs_step - 1) // abs_step) shape += (slice_size,) elif isinstance(i, (Sequence, np.ndarray)): if is_unpacked_axis: @@ -585,7 +621,7 @@ def get_index_parameters(self, index): unpacked_index += (slice(None),) else: if isinstance(index, np.ndarray) and index.dtype == np.bool: - mask_packed = [False] * self.shape[axis] + mask_packed = [False] * packed_shape[axis] for j, value in enumerate(i): mask_packed[j // bit_width] |= True packed_index = mask_packed @@ -594,7 +630,7 @@ def get_index_parameters(self, index): else: # adjust indexing for this packed axis data = np.array([s // bit_width for s in i], dtype=self.dtype) - # remove duplicate entries, including for nested arrays + # remove duplicate entries, including nested arrays if data.ndim > 1: rows = [] for j, row_data in enumerate(data): @@ -608,9 +644,9 @@ def get_index_parameters(self, index): data = data[np.sort(unique_indices)] packed_index += (data,) - shape += (self.shape[axis],) + shape += (packed_shape[axis],) if axes_in_index == 1: - unpacked_index += (i,) + unpacked_index += ([s % bit_width for s in i],) elif i is np.newaxis: packed_index += (i,) unpacked_index += (i,) @@ -624,10 +660,10 @@ def get_index_parameters(self, index): return packed_index, unpacked_index, shape - # TODO: Should this be the default shape returned and a cast (i.e. .view(np.ndarray) is needed to get the packed shape? @property - def unpacked_shape(self): - return self.shape[:-1] + (self._axis_count,) + def shape(self): + # A cast to np.ndarray is needed to get the packed shape + return self.view(np.ndarray).shape[:-1] + (self._axis_count,) def get_unpacked_value(self, index): # Numpy indexing is handled primarily in https://github.com/numpy/numpy/blob/maintenance/1.26.x/numpy/core/src/multiarray/mapping.c#L1435 From b2f316cff7ef69ed9863f0171e202a6f52fcb117 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Fri, 10 Jan 2025 16:28:23 -0800 Subject: [PATCH 20/21] Add additional tests; Fix multiply/outer product broadcasting --- src/galois/_fields/_gf2.py | 28 ++++++----- tests/fields/test_bitpacked.py | 85 ++++++++++++++++++++++++++++++---- 2 files changed, 89 insertions(+), 24 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 04fc16d3e..388be622f 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -172,22 +172,20 @@ def __call__(self, ufunc, method, inputs, kwargs, meta): else: result_shape = np.broadcast_shapes(*(i.shape for i in inputs)) - if is_outer_product and len(inputs) == 2: - a = np.unpackbits(inputs[0]) - output = a[:, np.newaxis].view(np.ndarray) * inputs[1].view(np.ndarray) + if is_outer_product: + assert len(inputs) == 2 + # Unpack the first argument and propagate the bitpacked second argument + inputs = [np.unpackbits(x).view(np.ndarray) if i == 0 else x.view(np.ndarray) for i, x in enumerate(inputs)] + output = np.multiply.outer(*inputs) else: - output = super().__call__(ufunc, method, inputs, kwargs, meta) - - assert len(output.view(np.ndarray).shape) == len(result_shape) - # output = output.view(np.ndarray) - # if output.shape != result_shape: - # for axis, shape in enumerate(zip(output.shape, result_shape)): - # if axis == len(result_shape) - 1: - # # The last axis remains packed - # break - # - # if shape[0] != shape[1]: - # output = np.unpackbits(output, axis=axis, count=shape[1]) + if any(i.shape != inputs[0].shape for i in inputs): + # We can't do simple bitwise multiplication when the shapes aren't the same due to broadcasting + inputs = [np.unpackbits(i) for i in inputs] + output = reduce(operator.mul, inputs) # We need this to use GF2's multiply + output = np.packbits(output) + else: + output = super().__call__(ufunc, method, inputs, kwargs, meta) + output = self.field._view(output) output._axis_count = result_shape[-1] return output diff --git a/tests/fields/test_bitpacked.py b/tests/fields/test_bitpacked.py index d304bed11..57a66e68c 100644 --- a/tests/fields/test_bitpacked.py +++ b/tests/fields/test_bitpacked.py @@ -1,5 +1,7 @@ import numpy as np import galois +from galois import GF2 + def test_galois_array_indexing(): # Define a Galois field array @@ -60,9 +62,6 @@ def test_galois_array_indexing(): # taken = np.take(arr, [0, 2]) # assert np.array_equal(taken, GF([1, 1])) - print("All tests passed.") - - def test_galois_array_setting(): # Define a Galois field array GF = galois.GF(2) @@ -134,9 +133,77 @@ def test_galois_array_setting(): arr_2d[np.ix_(row_indices, col_indices)] = GF([[0, 0], [0, 0]]) assert np.array_equal(arr_2d, np.packbits(GF([[0, 0], [0, 0]]))) - print("All set-indexing tests passed.") - - -if __name__ == "__main__": - test_galois_array_indexing() - test_galois_array_setting() +def test_inv(): + N = 10 + u = GF2.Random((N, N), seed=2) + p = np.packbits(u) + # print(x.get_unpacked_slice(1)) + # index = np.index_exp[:,1:4:2] + # index = np.index_exp[[0,1], [0, 1]] + # print(a) + # print(a[index]) + # print(x.get_unpacked_slice(index)) + print(np.linalg.inv(u)) + print(np.unpackbits(np.linalg.inv(p))) + assert np.array_equal(np.linalg.inv(u), np.unpackbits(np.linalg.inv(p))) + +def test_arithmetic(): + size = (20, 10) + cm = np.random.randint(2, size=size, dtype=np.uint8) + cm2 = np.random.randint(2, size=size, dtype=np.uint8) + vec = np.random.randint(2, size=size[1], dtype=np.uint8) + + cm_GF2 = GF2(cm) + cm2_GF2 = GF2(cm2) + cm3_GF2 = GF2(cm2.T) + vec_GF2 = GF2(vec) + + cm_gf2bp = np.packbits(cm_GF2) + cm2_gf2bp = np.packbits(cm2_GF2) + cm3_gf2bp = np.packbits(cm2_GF2.T) + vec_gf2bp = np.packbits(vec_GF2) + + # Addition + assert np.array_equal(np.unpackbits(cm_gf2bp + cm2_gf2bp), cm_GF2 + cm2_GF2) + + # Multiplication + assert np.array_equal(np.unpackbits(cm_gf2bp * cm2_gf2bp), cm_GF2 * cm2_GF2) + + # Matrix-vector product + assert np.array_equal(np.unpackbits(cm_gf2bp @ vec_gf2bp), cm_GF2 @ vec_GF2) + + # Matrix-matrix product + assert np.array_equal(np.unpackbits(cm_gf2bp @ cm3_gf2bp), cm_GF2 @ cm3_GF2) + +def test_broadcasting(): + a = np.random.randint(0, 2, 10) + b = np.random.randint(0, 2, 10) + x = np.packbits(GF2(a)) + y = np.packbits(GF2(b)) + + c = a * b + z = x * y + assert c.shape == z.shape == np.unpackbits(z).shape # (10,) + + c = np.multiply.outer(a, b) + z = np.multiply.outer(x, y) + assert np.array_equal(np.unpackbits(z), c) + assert c.shape == z.shape == np.unpackbits(z).shape # (10, 10) + +def test_advanced_broadcasting(): + a = np.random.randint(0, 2, (1, 2, 3)) + b = np.random.randint(0, 2, (2, 2, 1)) + x = np.packbits(GF2(a)) + y = np.packbits(GF2(b)) + + c = a * b + z = x * y + assert np.array_equal(np.unpackbits(z), c) + assert c.shape == z.shape == np.unpackbits(z).shape # (2, 2, 3) + + c = np.multiply.outer(a, b) + z = np.multiply.outer(x, y) + print(c.shape) + print(z.shape) + assert np.array_equal(np.unpackbits(z), c) + assert c.shape == z.shape == np.unpackbits(z).shape # (1, 2, 3, 2, 2, 1) From c535cd6d4f5fd52c864f08d346cc61054cafe7d4 Mon Sep 17 00:00:00 2001 From: Amir Ebrahimi Date: Fri, 10 Jan 2025 16:50:24 -0800 Subject: [PATCH 21/21] Fix up broadcasting in add, sub, div --- src/galois/_fields/_gf2.py | 33 ++++++++++++++++++++++++++------ tests/fields/test_bitpacked.py | 35 ++++++++++++++++++---------------- 2 files changed, 46 insertions(+), 22 deletions(-) diff --git a/src/galois/_fields/_gf2.py b/src/galois/_fields/_gf2.py index 388be622f..75b9a28c9 100644 --- a/src/galois/_fields/_gf2.py +++ b/src/galois/_fields/_gf2.py @@ -144,8 +144,15 @@ class add_ufunc_bitpacked(add_ufunc): """ def __call__(self, ufunc, method, inputs, kwargs, meta): - output = super().__call__(ufunc, method, inputs, kwargs, meta) - output._axis_count = inputs[0]._axis_count + result_shape = np.broadcast_shapes(*(i.shape for i in inputs)) + if any(i.shape != inputs[0].shape for i in inputs): + # We can't do simple bitwise addition when the shapes aren't the same due to broadcasting + inputs = [np.unpackbits(i) for i in inputs] + output = reduce(operator.add, inputs) # We need this to use GF2's addition + output = np.packbits(output) + else: + output = super().__call__(ufunc, method, inputs, kwargs, meta) + output._axis_count = result_shape[-1] return output @@ -155,8 +162,15 @@ class subtract_ufunc_bitpacked(subtract_ufunc): """ def __call__(self, ufunc, method, inputs, kwargs, meta): - output = super().__call__(ufunc, method, inputs, kwargs, meta) - output._axis_count = max(i._axis_count for i in inputs) + result_shape = np.broadcast_shapes(*(i.shape for i in inputs)) + if any(i.shape != inputs[0].shape for i in inputs): + # We can't do simple bitwise subtraction when the shapes aren't the same due to broadcasting + inputs = [np.unpackbits(i) for i in inputs] + output = reduce(operator.sub, inputs) # We need this to use GF2's subtraction + output = np.packbits(output) + else: + output = super().__call__(ufunc, method, inputs, kwargs, meta) + output._axis_count = result_shape[-1] return output @@ -197,8 +211,15 @@ class divide_ufunc_bitpacked(divide): """ def __call__(self, ufunc, method, inputs, kwargs, meta): - output = super().__call__(ufunc, method, inputs, kwargs, meta) - output._axis_count = max(i._axis_count for i in inputs) + result_shape = np.broadcast_shapes(*(i.shape for i in inputs)) + if any(i.shape != inputs[0].shape for i in inputs): + # We can't do simple bitwise division when the shapes aren't the same due to broadcasting + inputs = [np.unpackbits(i) for i in inputs] + output = reduce(operator.truediv, inputs) # We need this to use GF2's division + output = np.packbits(output) + else: + output = super().__call__(ufunc, method, inputs, kwargs, meta) + output._axis_count = result_shape[-1] return output diff --git a/tests/fields/test_bitpacked.py b/tests/fields/test_bitpacked.py index 57a66e68c..a8ca93ab2 100644 --- a/tests/fields/test_bitpacked.py +++ b/tests/fields/test_bitpacked.py @@ -1,6 +1,7 @@ import numpy as np import galois from galois import GF2 +import operator as ops def test_galois_array_indexing(): @@ -176,14 +177,15 @@ def test_arithmetic(): assert np.array_equal(np.unpackbits(cm_gf2bp @ cm3_gf2bp), cm_GF2 @ cm3_GF2) def test_broadcasting(): - a = np.random.randint(0, 2, 10) - b = np.random.randint(0, 2, 10) - x = np.packbits(GF2(a)) - y = np.packbits(GF2(b)) + a = GF2(np.random.randint(0, 2, 10)) + b = GF2(np.random.randint(0, 2, 10)) + x = np.packbits(a) + y = np.packbits(b) - c = a * b - z = x * y - assert c.shape == z.shape == np.unpackbits(z).shape # (10,) + for op in [ops.add, ops.sub, ops.mul]: + c = op(a, b) + z = op(x, y) + assert c.shape == z.shape == np.unpackbits(z).shape # (10,) c = np.multiply.outer(a, b) z = np.multiply.outer(x, y) @@ -191,15 +193,16 @@ def test_broadcasting(): assert c.shape == z.shape == np.unpackbits(z).shape # (10, 10) def test_advanced_broadcasting(): - a = np.random.randint(0, 2, (1, 2, 3)) - b = np.random.randint(0, 2, (2, 2, 1)) - x = np.packbits(GF2(a)) - y = np.packbits(GF2(b)) - - c = a * b - z = x * y - assert np.array_equal(np.unpackbits(z), c) - assert c.shape == z.shape == np.unpackbits(z).shape # (2, 2, 3) + a = GF2(np.random.randint(0, 2, (1, 2, 3))) + b = GF2(np.random.randint(0, 2, (2, 2, 1))) + x = np.packbits(a) + y = np.packbits(b) + + for op in [ops.add, ops.sub, ops.mul]: + c = op(a, b) + z = op(x, y) + assert np.array_equal(np.unpackbits(z), c) + assert c.shape == z.shape == np.unpackbits(z).shape # (2, 2, 3) c = np.multiply.outer(a, b) z = np.multiply.outer(x, y)