From ada4f11c6efd451a97265ae420cffe983042901a Mon Sep 17 00:00:00 2001 From: Erik Welch Date: Tue, 3 Nov 2026 12:33:48 -0800 Subject: [PATCH 1/9] Select by a transposed comparison as by its mask, not by the untransposed one select._match_expr rewrote A.select(A.T == x) as select.valueeq(A, x): its clause "args[0] is parent or args[0].T is parent" matched the transpose of the parent, and the rewrite then compared A itself, dropping the transpose and, for a matrix that is not square, the dimension check: S = Matrix.from_coo([0, 1], [1, 0], [1, 5], nrows=2, ncols=2) S.select(S.T == 1) before: (0, 1) = 1 after: (1, 0) = 5, as S.select((S.T == 1).new()) gives B (2 x 3): B.select(B.T == 1) before: a 2 x 3 result after: DimensionMismatch The clause is now "args[0] is parent"; anything else takes the mask path. On main since #252. --- graphblas/select/__init__.py | 2 +- graphblas/tests/test_matrix.py | 18 ++++++++++++++++++ 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/graphblas/select/__init__.py b/graphblas/select/__init__.py index df598667e..a2d94a0a0 100644 --- a/graphblas/select/__init__.py +++ b/graphblas/select/__init__.py @@ -116,7 +116,7 @@ def _match_expr(parent, expr): and len(args) == 2 and op.parent in _binary_to_select and args[1]._is_scalar - and (args[0] is parent or getattr(args[0], "T", None) is parent) + and args[0] is parent and op.type in _binary_to_select[op.parent] ): selectop = _binary_to_select[op.parent][op.type] diff --git a/graphblas/tests/test_matrix.py b/graphblas/tests/test_matrix.py index 2079897ec..f0fe1585a 100644 --- a/graphblas/tests/test_matrix.py +++ b/graphblas/tests/test_matrix.py @@ -1269,6 +1269,24 @@ def test_apply_indexunary(A): assert pickle.loads(pickle.dumps(indexunary.tril[int])) is indexunary.tril[int] +@autocompute +def test_select_by_a_transposed_comparison(): + """``A.select(A.T == x)`` keeps the entries where ``A.T == x``, as a computed mask does. + + It was rewritten to ``select.valueeq(A, x)``, which drops the transpose, and + for a matrix that is not square the DimensionMismatch (on main since #252). + """ + S = Matrix.from_coo([0, 1], [1, 0], [1, 5], nrows=2, ncols=2) + expected = Matrix.from_coo([1], [0], [5], nrows=2, ncols=2) + assert S.select(S.T == 1).new().isequal(expected) + assert S.select((S.T == 1).new()).new().isequal(expected) + # The rewrite still serves the comparison of the object itself + assert S.select(S == 1).new().isequal(Matrix.from_coo([0], [1], [1], nrows=2, ncols=2)) + B = Matrix.from_coo([0, 1], [1, 2], [1, 1], nrows=2, ncols=3) + with pytest.raises(DimensionMismatch): + B.select(B.T == 1).new() + + def test_select(A): A3 = Matrix.from_coo([0, 3, 3, 6], [3, 0, 2, 4], [3, 3, 3, 3], nrows=7, ncols=7) w1 = A.select(select.valueeq, 3).new() From 4a8afd2e4838a9e8de3d95d261122066edad39a4 Mon Sep 17 00:00:00 2001 From: Erik Welch Date: Wed, 4 Nov 2026 11:00:54 -0800 Subject: [PATCH 2/9] Clear an exported object's handle as soon as GraphBLAS frees it GraphBLAS frees the Matrix or Vector as it hands the arrays over in an export. _export gave GraphBLAS a copy of the handle and cleared the object's own handle only at the end, so an error in between, such as psg's claim_buffer_2d raising for an array UDT, left that handle pointing at the freed object: A.ss.export(format, give_ownership=True) raising after the export, then A.nvals before: UninitializedObject (a read of freed memory) after: NullPointer, as after an export that succeeds A plain export, of a copy, also freed that copy a second time when it was collected, which can abort the process (it did on SuiteSparse 7.4). GraphBLAS now gets the object's own handle, which it sets to NULL. --- graphblas/core/ss/matrix.py | 6 +++--- graphblas/core/ss/vector.py | 6 +++--- graphblas/tests/test_matrix.py | 29 +++++++++++++++++++++++++++++ 3 files changed, 35 insertions(+), 6 deletions(-) diff --git a/graphblas/core/ss/matrix.py b/graphblas/core/ss/matrix.py index 443019557..c68b1c626 100644 --- a/graphblas/core/ss/matrix.py +++ b/graphblas/core/ss/matrix.py @@ -863,7 +863,9 @@ def _export(self, format=None, *, sort=False, give_ownership=False, raw=False, m return rv if method == "export": - mhandle = ffi_new("GrB_Matrix*", parent._carg) + # GraphBLAS frees the Matrix as it hands the arrays over, and sets this + # handle to NULL, so an error after that cannot leave it dangling + mhandle = parent.gb_obj type_ = ffi_new("GrB_Type*") nrows_ = ffi_new("GrB_Index*") ncols_ = ffi_new("GrB_Index*") @@ -1170,8 +1172,6 @@ def _export(self, format=None, *, sort=False, give_ownership=False, raw=False, m rv["is_iso"] = is_iso rv["format"] = format rv["values"] = values - if method == "export": - parent.gb_obj[0] = NULL if parent.dtype._is_udt: rv["dtype"] = parent.dtype return rv diff --git a/graphblas/core/ss/vector.py b/graphblas/core/ss/vector.py index 8ee259510..4ad1c3eb9 100644 --- a/graphblas/core/ss/vector.py +++ b/graphblas/core/ss/vector.py @@ -515,7 +515,9 @@ def _export(self, format=None, *, sort=False, give_ownership=False, raw=False, m size = parent._size if method == "export": - vhandle = ffi_new("GrB_Vector*", parent._carg) + # GraphBLAS frees the Vector as it hands the arrays over, and sets this + # handle to NULL, so an error after that cannot leave it dangling + vhandle = parent.gb_obj type_ = ffi_new("GrB_Type*") size_ = ffi_new("GrB_Index*") args = (type_, size_) @@ -642,8 +644,6 @@ def _export(self, format=None, *, sort=False, give_ownership=False, raw=False, m format=format, values=values, ) - if method == "export": - parent.gb_obj[0] = NULL if parent.dtype._is_udt: rv["dtype"] = parent.dtype return rv diff --git a/graphblas/tests/test_matrix.py b/graphblas/tests/test_matrix.py index f0fe1585a..d24e82972 100644 --- a/graphblas/tests/test_matrix.py +++ b/graphblas/tests/test_matrix.py @@ -21,6 +21,7 @@ InvalidObject, InvalidValue, NotImplementedException, + NullPointer, OutputNotEmpty, ) @@ -251,6 +252,34 @@ def test_nvals(A): assert A.nvals == 12 +@pytest.mark.skipif("not suitesparse") +def test_ss_export_failing_after_graphblas_frees_the_object(monkeypatch): + # GraphBLAS frees the object as it hands its arrays over. An error after that, + # here in claiming the values, left the handle pointing at the freed object, so a + # later call read freed memory and the object was freed again when collected. + import gc + + import graphblas.core.ss.matrix as ss_matrix + import graphblas.core.ss.vector as ss_vector + + def fail(*args): + raise RuntimeError("claim failed") + + monkeypatch.setattr(ss_matrix, "claim_buffer_2d", fail) + monkeypatch.setattr(ss_vector, "claim_buffer", fail) + A = Matrix.from_dense([[1.5, 2.5], [3.5, 4.5]]) + v = Vector.from_dense([1.5, 2.5]) + for obj, format, nvals in [(A, "fullr", 4), (v, "full", 2)]: + with pytest.raises(RuntimeError, match="claim failed"): + obj.ss.export(format) # exports a copy + gc.collect() # frees the copy + assert obj.nvals == nvals + with pytest.raises(RuntimeError, match="claim failed"): + obj.ss.export(format, give_ownership=True) + with pytest.raises(NullPointer): + obj._nvals + + def test_build(A): assert A.nvals == 12 A.clear() From 29d56f7ea1bf3ba9ab333066ce03182aa39098b1 Mon Sep 17 00:00:00 2001 From: Erik Welch Date: Wed, 4 Nov 2026 11:17:08 -0800 Subject: [PATCH 3/9] Say an infinity is close only to itself in isclose, as math.isclose does isclose used math.isclose's formula without its infinity check. With an infinite x or y, abs(x - y) is inf and so is rel_tol * max(abs(x), abs(y)), so infinity was close to every number, in binary.isclose (Vector and Matrix isclose with Numba, Scalar.isclose) and in _isclose_recipe (without Numba): Scalar.from_value(inf).isclose(1.0), .isclose(-inf) Vector.from_coo([0], [inf]).isclose(Vector.from_coo([0], [1.0])) before: True (all three) after: False, as math.isclose and numpy.isclose say Two infinities of the same sign are equal, and so still close; NaN is still close to nothing. --- graphblas/core/matrix.py | 3 ++- graphblas/core/operator/binary.py | 13 +++++++++++-- graphblas/core/operator/unary.py | 6 ++++-- graphblas/core/scalar.py | 6 ++++-- graphblas/core/vector.py | 13 +++++++++++-- graphblas/tests/test_matrix.py | 17 +++++++++++++++++ graphblas/tests/test_scalar.py | 14 ++++++++++++++ graphblas/tests/test_vector.py | 17 +++++++++++++++++ 8 files changed, 80 insertions(+), 9 deletions(-) diff --git a/graphblas/core/matrix.py b/graphblas/core/matrix.py index f92a4b028..299d98d8c 100644 --- a/graphblas/core/matrix.py +++ b/graphblas/core/matrix.py @@ -459,7 +459,8 @@ def isequal(self, other, *, check_dtype=False, **opts): def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False, **opts): """Check for approximate equality (including same size and same structure). - Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``. + Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``, + except that, as in ``math.isclose``, an infinity is close only to itself. Parameters ---------- diff --git a/graphblas/core/operator/binary.py b/graphblas/core/operator/binary.py index 0db22748c..bb284c1a4 100644 --- a/graphblas/core/operator/binary.py +++ b/graphblas/core/operator/binary.py @@ -339,7 +339,12 @@ def _rpow(x, y): def _isclose(rel_tol=1e-7, abs_tol=0.0): def inner(x, y): # pragma: no cover (numba) - return x == y or abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol) + if x == y: + return True + # As in math.isclose, an infinity is close only to itself + if abs(x) == np.inf or abs(y) == np.inf: + return False + return abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol) return inner @@ -847,7 +852,11 @@ def register_new( >>> def user_isclose(rel_tol=1e-7, abs_tol=0.0): >>> def inner(x, y): - >>> return x == y or abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol) + >>> if x == y: + >>> return True + >>> if abs(x) == np.inf or abs(y) == np.inf: + >>> return False + >>> return abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol) >>> return inner >>> gb.binary.register_new("user_isclose", user_isclose, parameterized=True) diff --git a/graphblas/core/operator/unary.py b/graphblas/core/operator/unary.py index c0bdf0434..e8f0870e0 100644 --- a/graphblas/core/operator/unary.py +++ b/graphblas/core/operator/unary.py @@ -319,7 +319,8 @@ def register_anonymous(cls, func, name=None, *, parameterized=False, is_udt=Fals can be used as: ``A.ewise_mult(B, gb.binary.isclose(rel_tol=1e-5))``. When creating a parameterized user-defined operator, the ``func`` parameter must be a callable that *returns* a function that will - then get compiled. See the ``user_isclose`` example below. + then get compiled. See the ``user_isclose`` example in + ``BinaryOp.register_new``. is_udt : bool, default False Whether the operator is intended to operate on user-defined types. If True, then the function will not be automatically compiled for @@ -357,7 +358,8 @@ def register_new(cls, name, func, *, parameterized=False, is_udt=False, lazy=Fal can be used as: ``A.ewise_mult(B, gb.binary.isclose(rel_tol=1e-5))``. When creating a parameterized user-defined operator, the ``func`` parameter must be a callable that *returns* a function that will - then get compiled. See the ``user_isclose`` example below. + then get compiled. See the ``user_isclose`` example in + ``BinaryOp.register_new``. is_udt : bool, default False Whether the operator is intended to operate on user-defined types. If True, then the function will not be automatically compiled for diff --git a/graphblas/core/scalar.py b/graphblas/core/scalar.py index 5f1ac7c01..de7a2353c 100644 --- a/graphblas/core/scalar.py +++ b/graphblas/core/scalar.py @@ -17,6 +17,7 @@ _literal_type, ) from .operator import get_typed_op +from .operator.binary import _isclose from .udt import ( _element_bytes, _element_shapes_differ, @@ -34,7 +35,7 @@ if _supports_udfs: from ..binary import isclose else: - from .operator.binary import _isclose as isclose + isclose = _isclose ffi_new = ffi.new @@ -271,7 +272,8 @@ def isequal(self, other, *, check_dtype=False): def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False): """Check for approximate equality (including whether the value is missing). - Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``. + Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``, + except that, as in ``math.isclose``, an infinity is close only to itself. Parameters ---------- diff --git a/graphblas/core/vector.py b/graphblas/core/vector.py index ddb69fcdf..9e88fbeb0 100644 --- a/graphblas/core/vector.py +++ b/graphblas/core/vector.py @@ -111,7 +111,8 @@ def _select_mask(updater, obj, mask): def _isclose_recipe(self, other, rel_tol, abs_tol, **opts): - # x == y or abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol) + # x == y or abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol), except that + # an infinity is close only to itself isequal = self.ewise_mult(other, binary.eq).new(bool, name="isclose", **opts) if isequal._nvals != self._nvals: return False @@ -131,6 +132,13 @@ def _isclose_recipe(self, other, rel_tol, abs_tol, **opts): # max(abs(x), abs(y)) x(**opts) << x.ewise_mult(y, binary.max) max_x_y = x + # As in math.isclose, an infinity is close only to itself, and those are equal + if max_x_y.ndim == 1: + biggest = max_x_y.reduce(monoid.max, allow_empty=False).new(**opts).value + else: + biggest = max_x_y.reduce_scalar(monoid.max, allow_empty=False).new(**opts).value + if biggest == np.inf: + return False # rel_tol * max(abs(x), abs(y)) max_x_y(**opts) << max_x_y.apply(binary.times, rel_tol) # max(rel_tol * max(abs(x), abs(y)), abs_tol) @@ -398,7 +406,8 @@ def isequal(self, other, *, check_dtype=False, **opts): def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False, **opts): """Check for approximate equality (including same size and same structure). - Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``. + Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``, + except that, as in ``math.isclose``, an infinity is close only to itself. Parameters ---------- diff --git a/graphblas/tests/test_matrix.py b/graphblas/tests/test_matrix.py index d24e82972..238941eef 100644 --- a/graphblas/tests/test_matrix.py +++ b/graphblas/tests/test_matrix.py @@ -1,5 +1,6 @@ import inspect import itertools +import math import pickle import sys import types @@ -1982,6 +1983,22 @@ def test_isclose(A, v): assert C6.isclose(A, rel_tol=1e-3) +def test_isclose_infinity(): + # As in math.isclose, an infinity is close only to itself. It was close to every + # number (abs(inf - y) is inf, and so is rel_tol * inf), with Numba and without. + from graphblas.core.vector import _isclose_recipe + + values = [math.inf, -math.inf, 1.0, 0.0, math.nan] + for x, y in itertools.product(values, repeat=2): + left = Matrix.from_coo([0, 1], [1, 0], [x, 5.0]) + right = Matrix.from_coo([0, 1], [1, 0], [y, 5.0]) + for abs_tol in [0.0, math.inf]: + expected = math.isclose(x, y, rel_tol=1e-7, abs_tol=abs_tol) + if supports_udfs: + assert left.isclose(right, abs_tol=abs_tol) is expected, (x, y) + assert _isclose_recipe(left, right, 1e-7, abs_tol) is expected, (x, y) + + @pytest.mark.slow def test_transpose_equals(A): data = [ diff --git a/graphblas/tests/test_scalar.py b/graphblas/tests/test_scalar.py index 2d5974032..5bf6408f8 100644 --- a/graphblas/tests/test_scalar.py +++ b/graphblas/tests/test_scalar.py @@ -1,4 +1,6 @@ import inspect +import itertools +import math import pickle import random import sys @@ -199,6 +201,18 @@ def test_isclose(): assert Scalar.from_value(None, dtype="FP64").isclose(Scalar.from_value(None, dtype="FP32")) +def test_isclose_infinity(): + # As in math.isclose, an infinity is close only to itself. It was close to every + # number (abs(inf - y) is inf, and so is rel_tol * inf), in both of isclose's paths. + values = [math.inf, -math.inf, 1.0, 0.0, math.nan] + for x, y in itertools.product(values, repeat=2): + for abs_tol in [0.0, math.inf]: + expected = math.isclose(x, y, rel_tol=1e-7, abs_tol=abs_tol) + assert Scalar.from_value(x).isclose(y, abs_tol=abs_tol) is expected, (x, y) + other = Scalar.from_value(y) + assert Scalar.from_value(x).isclose(other, abs_tol=abs_tol) is expected, (x, y) + + def test_isclose_udt_raises(): """Reject ``isclose`` on UDTs with a clear message that points to isequal.""" arr_udt = dtypes.register_anonymous(np.dtype((np.int64, (3,))), "_IscloseUdt3") diff --git a/graphblas/tests/test_vector.py b/graphblas/tests/test_vector.py index 541a892f9..d25bce0a7 100644 --- a/graphblas/tests/test_vector.py +++ b/graphblas/tests/test_vector.py @@ -1,5 +1,6 @@ import inspect import itertools +import math import pickle import sys import types @@ -1137,6 +1138,22 @@ def test_isclose(v): assert not u4.isclose(u4b) +def test_isclose_infinity(): + # As in math.isclose, an infinity is close only to itself. It was close to every + # number (abs(inf - y) is inf, and so is rel_tol * inf), with Numba and without. + from graphblas.core.vector import _isclose_recipe + + values = [math.inf, -math.inf, 1.0, 0.0, math.nan] + for x, y in itertools.product(values, repeat=2): + left = Vector.from_coo([0, 2], [x, 5.0], size=3) + right = Vector.from_coo([0, 2], [y, 5.0], size=3) + for abs_tol in [0.0, math.inf]: + expected = math.isclose(x, y, rel_tol=1e-7, abs_tol=abs_tol) + if supports_udfs: + assert left.isclose(right, abs_tol=abs_tol) is expected, (x, y) + assert _isclose_recipe(left, right, 1e-7, abs_tol) is expected, (x, y) + + def test_binary_op(v): v2 = v.dup() v2[1] = 0 From da33b3fa31c2f18683cdba175e2892fbadf27c9b Mon Sep 17 00:00:00 2001 From: Erik Welch Date: Sat, 7 Nov 2026 20:54:04 -0600 Subject: [PATCH 4/9] Compute isclose on integers without wrapping in x - y or abs(x) isclose computed abs(x - y) and abs(x) in the integer type, where they wrap: in UINT64, and with Numba in uint64 for any unsigned type (Numba widens to 64 bits); in INT64 at its extremes, as is abs(-2**63); and without Numba in every integer type. So, for Vector, Matrix and Scalar: UINT64 [0] and [2**64 - 1]: before True, after False UINT8 or UINT32 5 and 6, rel_tol=0.5: before False, after True INT64 -2**63 and 2**63 - 1: before True, after False INT8 -128 and 127, rel_tol=0.5: before True without Numba, after False and Scalar.isclose raised OverflowError for a UINT64 value above INT64's maximum, because Numba typed the Python int as int64. The difference is now exact (in uint64, for an integer type), the magnitudes are in FP64, and Scalar.isclose passes its values in the op's input types, as ewise_mult does. Mixed types are compared in the type they unify to, as before, so INT64 with UINT64 is in FP64. --- graphblas/core/operator/binary.py | 40 +++++++++++++++++++++++++++++-- graphblas/core/scalar.py | 6 ++++- graphblas/core/vector.py | 27 ++++++++++++++------- graphblas/tests/test_matrix.py | 16 +++++++++++++ graphblas/tests/test_scalar.py | 13 ++++++++++ graphblas/tests/test_vector.py | 31 ++++++++++++++++++++++++ 6 files changed, 121 insertions(+), 12 deletions(-) diff --git a/graphblas/core/operator/binary.py b/graphblas/core/operator/binary.py index bb284c1a4..6abccb39e 100644 --- a/graphblas/core/operator/binary.py +++ b/graphblas/core/operator/binary.py @@ -344,11 +344,46 @@ def inner(x, y): # pragma: no cover (numba) # As in math.isclose, an infinity is close only to itself if abs(x) == np.inf or abs(y) == np.inf: return False - return abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol) + return _abs_diff(x, y) <= max(rel_tol * max(_magnitude(x), _magnitude(y)), abs_tol) return inner +# Integer arithmetic in Numba wraps: x - y in int64 or uint64 (Numba widens smaller +# unsigned types to uint64, so 5 - 6 is 2**64 - 1), and abs at a signed type's +# minimum. Python ints do not, so for them these are plain Python; the overloads +# below compile them for Numba's integers so that they do not wrap either. +def _abs_diff(x, y): + return abs(x - y) + + +def _magnitude(x): + return abs(x) + + +if _has_numba: + from numba.extending import overload + + @overload(_abs_diff) + def _abs_diff_overload(x, y): + if isinstance(x, numba.types.Integer) and x == y: + + def abs_diff(x, y): # pragma: no cover (numba) + # The difference fits uint64, where both convert modulo 2**64 + if x >= y: + return np.uint64(x) - np.uint64(y) + return np.uint64(y) - np.uint64(x) + + return abs_diff + return lambda x, y: abs(x - y) + + @overload(_magnitude) + def _magnitude_overload(x): + if isinstance(x, numba.types.Integer): + return lambda x: abs(np.float64(x)) + return lambda x: abs(x) # noqa: PLW0108 (an overload returns a Python function) + + _MAX_INT64 = np.iinfo(np.int64).max @@ -848,7 +883,8 @@ def register_new( >>> dir(gb.binary) [..., 'max_zero', ...] - This is how ``gb.binary.isclose`` is defined: + ``gb.binary.isclose`` is defined like this, except that it computes ``x - y`` + and ``abs(x)`` for integers without wrapping: >>> def user_isclose(rel_tol=1e-7, abs_tol=0.0): >>> def inner(x, y): diff --git a/graphblas/core/scalar.py b/graphblas/core/scalar.py index de7a2353c..b3b78193d 100644 --- a/graphblas/core/scalar.py +++ b/graphblas/core/scalar.py @@ -338,7 +338,11 @@ def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False): is_right_scalar=True, kind="binary", ) - return isclose_func._numba_func(self.value, other.value) + # Pass the values in the op's input types, as ewise_mult does; Numba types a + # Python int as int64, which a UINT64 value need not fit + x = isclose_func.type.np_type.type(self.value) + y = isclose_func.type2.np_type.type(other.value) + return isclose_func._numba_func(x, y) def clear(self): """In-place operation which clears the value in the Scalar. diff --git a/graphblas/core/vector.py b/graphblas/core/vector.py index 9e88fbeb0..6e206f1b8 100644 --- a/graphblas/core/vector.py +++ b/graphblas/core/vector.py @@ -3,7 +3,7 @@ import numpy as np from .. import backend, binary, monoid, select, semiring, unary -from ..dtypes import _INDEX, FP64, INT64, lookup_dtype, unify +from ..dtypes import _INDEX, FP64, INT64, UINT64, lookup_dtype, unify from ..exceptions import DimensionMismatch, NoValue, check_status from . import _supports_udfs, automethods, ffi, lib, utils from .base import BaseExpression, BaseType, call @@ -110,6 +110,10 @@ def _select_mask(updater, obj, mask): updater << obj.dup(mask=mask) +def _abs_type(dtype): + return dtype if dtype.np_type.kind == "c" else FP64 + + def _isclose_recipe(self, other, rel_tol, abs_tol, **opts): # x == y or abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol), except that # an infinity is close only to itself @@ -125,10 +129,9 @@ def _isclose_recipe(self, other, rel_tol, abs_tol, **opts): # So we can use structural mask below isequal(**opts) << select.value(isequal == True) # noqa: E712 - # abs(x) - x = self.apply(unary.abs).new(FP64, mask=~isequal.S, **opts) - # abs(y) - y = other.apply(unary.abs).new(FP64, mask=~isequal.S, **opts) + # abs(x) and abs(y), in FP64 for a real type: abs wraps at an integer type's minimum + x = self.apply(unary.abs[_abs_type(self.dtype)]).new(FP64, mask=~isequal.S, **opts) + y = other.apply(unary.abs[_abs_type(other.dtype)]).new(FP64, mask=~isequal.S, **opts) # max(abs(x), abs(y)) x(**opts) << x.ewise_mult(y, binary.max) max_x_y = x @@ -144,11 +147,17 @@ def _isclose_recipe(self, other, rel_tol, abs_tol, **opts): # max(rel_tol * max(abs(x), abs(y)), abs_tol) max_x_y(**opts) << max_x_y.apply(binary.max, abs_tol) - # x - y - y(~isequal.S, replace=True, **opts) << self.ewise_mult(other, binary.minus) - abs_x_y = y # abs(x - y) - abs_x_y(**opts) << abs_x_y.apply(unary.abs) + if unify(self.dtype, other.dtype).np_type.kind in "iu": + # x - y wraps in an integer type, and the difference may not fit it. It fits + # UINT64, where max(x, y) and min(x, y) convert modulo 2**64. + hi = self.ewise_mult(other, binary.max).new(UINT64, mask=~isequal.S, **opts) + lo = self.ewise_mult(other, binary.min).new(UINT64, mask=~isequal.S, **opts) + y(~isequal.S, replace=True, **opts) << hi.ewise_mult(lo, binary.minus) + else: + y(~isequal.S, replace=True, **opts) << self.ewise_mult(other, binary.minus) + y(**opts) << y.apply(unary.abs) + abs_x_y = y # abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol) isequal(**opts) << abs_x_y.ewise_mult(max_x_y, binary.le) diff --git a/graphblas/tests/test_matrix.py b/graphblas/tests/test_matrix.py index 238941eef..8848de7dd 100644 --- a/graphblas/tests/test_matrix.py +++ b/graphblas/tests/test_matrix.py @@ -1999,6 +1999,22 @@ def test_isclose_infinity(): assert _isclose_recipe(left, right, 1e-7, abs_tol) is expected, (x, y) +def test_isclose_integers_do_not_wrap(): + # x - y or abs(x) wrapped in an integer type; see test_vector.ISCLOSE_WRAPPING_CASES + from graphblas.core.vector import _isclose_recipe + + from .test_vector import ISCLOSE_WRAPPING_CASES + + for xtype, ytype, x, y, rel_tol in ISCLOSE_WRAPPING_CASES: + expected = math.isclose(x, y, rel_tol=rel_tol) + for (lt, lv), (rt, rv) in [((xtype, x), (ytype, y)), ((ytype, y), (xtype, x))]: + left = Matrix.from_coo([0, 1], [1, 0], [lv, 5], dtype=lt) + right = Matrix.from_coo([1, 0], [0, 1], [rv, 5], dtype=rt).T + if supports_udfs: + assert left.isclose(right, rel_tol=rel_tol) is expected, (lt, lv, rt, rv) + assert _isclose_recipe(left, right, rel_tol, 0.0) is expected, (lt, lv, rt, rv) + + @pytest.mark.slow def test_transpose_equals(A): data = [ diff --git a/graphblas/tests/test_scalar.py b/graphblas/tests/test_scalar.py index 5bf6408f8..d4168cd21 100644 --- a/graphblas/tests/test_scalar.py +++ b/graphblas/tests/test_scalar.py @@ -213,6 +213,19 @@ def test_isclose_infinity(): assert Scalar.from_value(x).isclose(other, abs_tol=abs_tol) is expected, (x, y) +def test_isclose_integers_do_not_wrap(): + # x - y or abs(x) wrapped in an integer type, and a UINT64 value above INT64's + # maximum raised OverflowError with Numba; see test_vector.ISCLOSE_WRAPPING_CASES + from .test_vector import ISCLOSE_WRAPPING_CASES + + for xtype, ytype, x, y, rel_tol in ISCLOSE_WRAPPING_CASES: + expected = math.isclose(x, y, rel_tol=rel_tol) + for (lt, lv), (rt, rv) in [((xtype, x), (ytype, y)), ((ytype, y), (xtype, x))]: + left = Scalar.from_value(lv, dtype=lt) + right = Scalar.from_value(rv, dtype=rt) + assert left.isclose(right, rel_tol=rel_tol) is expected, (lt, lv, rt, rv) + + def test_isclose_udt_raises(): """Reject ``isclose`` on UDTs with a clear message that points to isequal.""" arr_udt = dtypes.register_anonymous(np.dtype((np.int64, (3,))), "_IscloseUdt3") diff --git a/graphblas/tests/test_vector.py b/graphblas/tests/test_vector.py index d25bce0a7..a1e410076 100644 --- a/graphblas/tests/test_vector.py +++ b/graphblas/tests/test_vector.py @@ -1154,6 +1154,37 @@ def test_isclose_infinity(): assert _isclose_recipe(left, right, 1e-7, abs_tol) is expected, (x, y) +# (x's dtype, y's dtype, x, y, rel_tol): each was wrong in at least one of isclose's +# paths, because x - y or abs(x) wrapped in an integer type. The values are far enough +# apart (or exactly on the tolerance) that FP64 gives math.isclose's answer too. +ISCLOSE_WRAPPING_CASES = [ + ("UINT64", "UINT64", 0, 2**64 - 1, 1e-7), + ("UINT64", "UINT64", 2**63, 2**64 - 1, 1e-7), + ("UINT32", "UINT32", 5, 6, 0.5), + ("UINT8", "UINT8", 5, 6, 0.5), + ("UINT8", "UINT64", 0, 2**64 - 1, 1e-7), + ("INT64", "INT64", -(2**63), 2**63 - 1, 1e-7), + ("INT64", "INT64", -(2**63), 0, 1.0), + ("INT8", "INT8", -128, 127, 0.5), + ("INT8", "INT8", -128, 0, 1.0), + ("INT8", "INT64", 0, -(2**63), 1e-7), + ("INT8", "UINT64", -1, 2**64 - 1, 1e-7), +] + + +def test_isclose_integers_do_not_wrap(): + from graphblas.core.vector import _isclose_recipe + + for xtype, ytype, x, y, rel_tol in ISCLOSE_WRAPPING_CASES: + expected = math.isclose(x, y, rel_tol=rel_tol) + for (lt, lv), (rt, rv) in [((xtype, x), (ytype, y)), ((ytype, y), (xtype, x))]: + left = Vector.from_coo([0, 2], [lv, 5], dtype=lt, size=3) + right = Vector.from_coo([0, 2], [rv, 5], dtype=rt, size=3) + if supports_udfs: + assert left.isclose(right, rel_tol=rel_tol) is expected, (lt, lv, rt, rv) + assert _isclose_recipe(left, right, rel_tol, 0.0) is expected, (lt, lv, rt, rv) + + def test_binary_op(v): v2 = v.dup() v2[1] = 0 From a2b0aecfb1074777dfb5ebdc93ed5b93eb579c60 Mon Sep 17 00:00:00 2001 From: Erik Welch Date: Sat, 7 Nov 2026 20:54:55 -0600 Subject: [PATCH 5/9] Keep the imaginary part of x - y in isclose without Numba _isclose_recipe, which Vector and Matrix isclose use without Numba, stored x - y in an FP64 Vector before taking abs, and so kept only the real part of a complex difference: Vector.from_coo([0], [0j]).isclose(Vector.from_coo([0], [1j])) before: True without Numba (False with it) after: False, as cmath.isclose says abs is now taken in the complex type first. abs(x) and abs(y) are taken in the type x - y is in, as the Numba op does, so an FC32 operand beside an FP64 one is not rounded to FP32 first. --- graphblas/core/vector.py | 13 ++++++++----- graphblas/tests/test_vector.py | 25 +++++++++++++++++++++++++ 2 files changed, 33 insertions(+), 5 deletions(-) diff --git a/graphblas/core/vector.py b/graphblas/core/vector.py index 6e206f1b8..82fb5270f 100644 --- a/graphblas/core/vector.py +++ b/graphblas/core/vector.py @@ -129,9 +129,11 @@ def _isclose_recipe(self, other, rel_tol, abs_tol, **opts): # So we can use structural mask below isequal(**opts) << select.value(isequal == True) # noqa: E712 - # abs(x) and abs(y), in FP64 for a real type: abs wraps at an integer type's minimum - x = self.apply(unary.abs[_abs_type(self.dtype)]).new(FP64, mask=~isequal.S, **opts) - y = other.apply(unary.abs[_abs_type(other.dtype)]).new(FP64, mask=~isequal.S, **opts) + # abs(x) and abs(y), in the type x - y is in below, or in FP64 for a real type: + # abs wraps at an integer type's minimum + abs_type = _abs_type(unify(self.dtype, other.dtype)) + x = self.apply(unary.abs[abs_type]).new(FP64, mask=~isequal.S, **opts) + y = other.apply(unary.abs[abs_type]).new(FP64, mask=~isequal.S, **opts) # max(abs(x), abs(y)) x(**opts) << x.ewise_mult(y, binary.max) max_x_y = x @@ -155,8 +157,9 @@ def _isclose_recipe(self, other, rel_tol, abs_tol, **opts): lo = self.ewise_mult(other, binary.min).new(UINT64, mask=~isequal.S, **opts) y(~isequal.S, replace=True, **opts) << hi.ewise_mult(lo, binary.minus) else: - y(~isequal.S, replace=True, **opts) << self.ewise_mult(other, binary.minus) - y(**opts) << y.apply(unary.abs) + # abs before storing in y, whose FP64 would keep only a complex number's real part + diff = self.ewise_mult(other, binary.minus).new(mask=~isequal.S, **opts) + y(~isequal.S, replace=True, **opts) << diff.apply(unary.abs) abs_x_y = y # abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol) diff --git a/graphblas/tests/test_vector.py b/graphblas/tests/test_vector.py index a1e410076..509b74173 100644 --- a/graphblas/tests/test_vector.py +++ b/graphblas/tests/test_vector.py @@ -1,3 +1,4 @@ +import cmath import inspect import itertools import math @@ -1185,6 +1186,30 @@ def test_isclose_integers_do_not_wrap(): assert _isclose_recipe(left, right, rel_tol, 0.0) is expected, (lt, lv, rt, rv) +@pytest.mark.skipif("not dtypes._supports_complex") +def test_isclose_complex(): + # Without Numba, isclose kept only the real part of x - y, so 0 was close to 1j + from graphblas.core.vector import _isclose_recipe + + values = [0j, 1j, -1j, 1 + 1j, 3 + 4j, 3 + 4.000001j] + for x, y in itertools.product(values, repeat=2): + expected = cmath.isclose(x, y, rel_tol=1e-6) + left = Vector.from_coo([0, 2], [x, 5], dtype="FC64", size=3) + right = Vector.from_coo([0, 2], [y, 5], dtype="FC64", size=3) + if supports_udfs: + assert left.isclose(right, rel_tol=1e-6) is expected, (x, y) + assert _isclose_recipe(left, right, 1e-6, 0.0) is expected, (x, y) + # abs(x) and abs(y) are in the type x - y is in, here FC64; abs(1 + 1j) in FC32 is + # smaller than in FC64, so the case exactly on the tolerance was not close + left = Vector.from_coo([0], [0.0]) + right = Vector.from_coo([0], [1 + 1j], dtype="FC32") + expected = cmath.isclose(0.0, complex(np.complex64(1 + 1j)), rel_tol=1.0) + assert expected + if supports_udfs: + assert left.isclose(right, rel_tol=1.0) is expected + assert _isclose_recipe(left, right, 1.0, 0.0) is expected + + def test_binary_op(v): v2 = v.dup() v2[1] = 0 From 634d4baa989574f6f5473eeae4e8abc5d5059705 Mon Sep 17 00:00:00 2001 From: Erik Welch Date: Wed, 4 Nov 2026 12:20:42 -0800 Subject: [PATCH 6/9] Count each row's entries when importing coo, so empty rows keep indptr whole import_coor and import_cooc build CSR (or CSC) index pointers with indices_to_indptr, which set indptr only at the rows that appear. Every empty row after the first one left a 0 in indptr, so an empty row in the middle moved entries to another row, and an empty last row lost every entry. Exporting as "coo" with sort=True, as "coor" or as "cooc" and importing back hit it, on main since 2021: A = Matrix.from_coo([0, 2], [1, 0], [1.0, 2.0], nrows=3, ncols=2) B = Matrix.ss.import_any(**A.ss.export("coo", sort=True)) before: B holds (2, 1) and (2, 0), and A.isequal(B) said True after: B holds (0, 1) and (2, 0), as A does A = Matrix.from_coo([0, 0], [0, 1], [1.0, 2.0], nrows=3, ncols=2) before: the round trip gives an empty Matrix after: both entries It now counts the entries of each row, then sums them. The indices are converted to uint64 first, so lists work, and one at or beyond nrows (or ncols) raises IndexOutOfBound, where it was written past the end of indptr with no bounds check under Numba. --- graphblas/core/ss/matrix.py | 37 +++++++++++++++++++++++++--------- graphblas/tests/test_matrix.py | 31 ++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 10 deletions(-) diff --git a/graphblas/core/ss/matrix.py b/graphblas/core/ss/matrix.py index c68b1c626..326afc21b 100644 --- a/graphblas/core/ss/matrix.py +++ b/graphblas/core/ss/matrix.py @@ -7,7 +7,13 @@ from ... import binary, monoid from ...dtypes import _INDEX, BOOL, INT64, UINT64, lookup_dtype -from ...exceptions import InvalidValue, _error_code_lookup, check_status, check_status_carg +from ...exceptions import ( + IndexOutOfBound, + InvalidValue, + _error_code_lookup, + check_status, + check_status_carg, +) from .. import NULL, _has_numba, ffi, lib from ..base import call from ..operator import get_typed_op @@ -3171,7 +3177,7 @@ def _import_coor( raise ValueError(f"Invalid format: {format!r}. Must be None or 'coor'.") if not sorted_rows: raise ValueError("sorted_rows must be True when importing 'coor' format") - indptr = indices_to_indptr(rows, nrows + 1) + indptr = _coo_indptr(rows, nrows, "rows", "nrows") return cls._import_csr( nrows=nrows, ncols=ncols, @@ -3344,7 +3350,7 @@ def _import_cooc( raise ValueError(f"Invalid format: {format!r}. Must be None or 'cooc'.") if not sorted_cols: raise ValueError("sorted_cols must be True when importing 'cooc' format") - indptr = indices_to_indptr(cols, ncols + 1) + indptr = _coo_indptr(cols, ncols, "cols", "ncols") return cls._import_csc( nrows=nrows, ncols=ncols, @@ -4558,19 +4564,30 @@ def issorted(arr): # pragma: no cover (numba) @njit def indices_to_indptr(indices, size): # pragma: no cover (numba) - """Calculate the indptr for e.g. CSR from sorted COO rows.""" + """Calculate the indptr for e.g. CSR from sorted uint64 COO rows, each below size - 1. + + Each row's count goes at ``indptr[row + 1]``, then a running sum. Setting only the + rows that appear left 0 for every empty row after the first one, so an empty row + in the middle moved entries to another row, and an empty last row lost them all. + """ indptr = np.zeros(size, dtype=indices.dtype) - index = np.uint64(0) one = np.uint64(1) for i in range(indices.size): - row = indices[i] - if row != index: - indptr[index + one] = i - index = row - indptr[index + one] = indices.size + indptr[indices[i] + one] += one + for i in range(1, size): + indptr[i] += indptr[i - 1] return indptr +def _coo_indptr(indices, size, name, size_name): + """Return the indptr, of ``size + 1``, of sorted COO ``indices`` (rows or columns).""" + indices = ints_to_numpy_buffer(indices, np.uint64, name=name) + # indices_to_indptr writes at each index + 1, with no bounds check under Numba + if indices.size > 0 and (biggest := indices.max()) >= size: + raise IndexOutOfBound(f"{name} must be less than {size_name}, {size}; got {biggest}") + return indices_to_indptr(indices, size + 1) + + @njit(parallel=True) def indptr_to_indices(indptr): # pragma: no cover (numba) indices = np.empty(indptr[-1], dtype=indptr.dtype) diff --git a/graphblas/tests/test_matrix.py b/graphblas/tests/test_matrix.py index 8848de7dd..8361030bd 100644 --- a/graphblas/tests/test_matrix.py +++ b/graphblas/tests/test_matrix.py @@ -4626,6 +4626,37 @@ def test_subarray_dtype_values_shape(): Matrix.from_coo([0, 1], [0, 1], np.ones((3, 2)), "FP64[2]") +@pytest.mark.skipif("not suitesparse") +def test_ss_import_coo_with_empty_rows_and_columns(): + # indptr was set only at the rows that appear, so every empty row after the first + # left a 0 in it: an empty row in the middle moved entries to another row (and + # isequal still said True), and an empty last row lost every entry + def entries(M): + return sorted(zip(*(x.tolist() for x in M.to_coo()), strict=True)) + + cases = [ + ([0, 2], [1, 0], 3, 2), # row 1, in the middle + ([0, 0], [0, 1], 3, 2), # rows 1 and 2, at the end + ([2, 2], [0, 1], 3, 2), # rows 0 and 1, at the start (this one worked) + ([1, 3], [2, 2], 5, 4), # rows and columns empty at the start, middle and end + ([], [], 3, 2), # no entries + ] + for rows, cols, nrows, ncols in cases: + A = Matrix.from_coo(rows, cols, np.arange(1.0, len(rows) + 1), nrows=nrows, ncols=ncols) + for M in [A, A.T.new()]: + for info in [M.ss.export("coor"), M.ss.export("cooc"), M.ss.export("coo", sort=True)]: + B = Matrix.ss.import_any(**info) + assert B.shape == M.shape + assert entries(B) == entries(M), (rows, cols, info["format"]) + # Lists work, and an index beyond the shape is refused before indptr is written + B = Matrix.ss.import_coor(rows=[0, 2], cols=[1, 0], values=[1.0, 2.0], nrows=3, ncols=2) + assert B.to_coo()[0].tolist() == [0, 2] + with pytest.raises(IndexOutOfBound, match="rows must be less than nrows, 3; got 3"): + Matrix.ss.import_coor(rows=[0, 3], cols=[0, 0], values=[1.0, 2.0], nrows=3, ncols=1) + with pytest.raises(IndexOutOfBound, match="cols must be less than ncols, 2; got 5"): + Matrix.ss.import_cooc(rows=[0, 0], cols=[1, 5], values=[1.0, 2.0], nrows=1, ncols=2) + + @pytest.mark.skipif("not suitesparse") def test_ss_import_values_size(): # GraphBLAS before 10.5.2 compares the byte size of imported values with the number From 72f79af6d8917a69b287a1c81c3af400c8e6590d Mon Sep 17 00:00:00 2001 From: Erik Welch Date: Sat, 7 Nov 2026 21:03:31 -0600 Subject: [PATCH 7/9] Read 2-D coo indices as a flat buffer, and count them with numpy without Numba The new indptr kernel for import_coor and import_cooc indexed a 2-D rows (or cols) array by its rows, past the end of its first axis: rows = np.array([[0, 1], [1, 2]], np.uint64) Matrix.ss.import_coor(rows=rows, cols=[0, 0, 1, 1], values=[1, 2, 3, 4], nrows=3, ncols=2) before: Segmentation fault with Numba; IndexError without it after: the Matrix from_coo(rows.ravel(), ...) gives, as other imports read 2-D indices Without Numba the kernel ran in Python, once per entry and once per row (1.1 s for 10**7 rows and 10**5 entries, where the old kernel took 12 ms). It now counts with np.bincount and sums with np.cumsum (20 ms). --- graphblas/core/ss/matrix.py | 14 ++++++++++++-- graphblas/tests/test_matrix.py | 7 +++++++ 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/graphblas/core/ss/matrix.py b/graphblas/core/ss/matrix.py index 326afc21b..3205a1a69 100644 --- a/graphblas/core/ss/matrix.py +++ b/graphblas/core/ss/matrix.py @@ -4581,11 +4581,21 @@ def indices_to_indptr(indices, size): # pragma: no cover (numba) def _coo_indptr(indices, size, name, size_name): """Return the indptr, of ``size + 1``, of sorted COO ``indices`` (rows or columns).""" - indices = ints_to_numpy_buffer(indices, np.uint64, name=name) + # Read 2-D indices as a flat buffer, as other imports do. The kernel would index + # a 2-D array by its rows, past the end of its first axis. + indices = ints_to_numpy_buffer(indices, np.uint64, name=name).ravel() # indices_to_indptr writes at each index + 1, with no bounds check under Numba if indices.size > 0 and (biggest := indices.max()) >= size: raise IndexOutOfBound(f"{name} must be less than {size_name}, {size}; got {biggest}") - return indices_to_indptr(indices, size + 1) + if _has_numba: + return indices_to_indptr(indices, size + 1) + # Without Numba the kernel's loops run in Python, once per entry and once per row. + # Every index is below size, so it reads the same as int64; counts are never negative. + indptr = np.empty(size + 1, np.uint64) + indptr[0] = 0 + counts = np.bincount(indices.view(np.int64), minlength=size) + np.cumsum(counts.view(np.uint64), out=indptr[1:]) + return indptr @njit(parallel=True) diff --git a/graphblas/tests/test_matrix.py b/graphblas/tests/test_matrix.py index 8361030bd..2f10e5ed0 100644 --- a/graphblas/tests/test_matrix.py +++ b/graphblas/tests/test_matrix.py @@ -4655,6 +4655,13 @@ def entries(M): Matrix.ss.import_coor(rows=[0, 3], cols=[0, 0], values=[1.0, 2.0], nrows=3, ncols=1) with pytest.raises(IndexOutOfBound, match="cols must be less than ncols, 2; got 5"): Matrix.ss.import_cooc(rows=[0, 0], cols=[1, 5], values=[1.0, 2.0], nrows=1, ncols=2) + # 2-D indices are a flat buffer, as in from_coo; they read past the end (or crashed) + rows = np.array([[0, 1], [1, 2]], np.uint64) + B = Matrix.ss.import_coor(rows=rows, cols=[0, 0, 1, 1], values=[1, 2, 3, 4], nrows=3, ncols=2) + assert entries(B) == entries(Matrix.from_coo(rows.ravel(), [0, 0, 1, 1], [1, 2, 3, 4])) + cols = np.array([[0], [0], [1]], np.uint64) + B = Matrix.ss.import_cooc(rows=[0, 2, 1], cols=cols, values=[1, 2, 3], nrows=3, ncols=2) + assert entries(B) == entries(Matrix.from_coo([0, 2, 1], cols.ravel(), [1, 2, 3])) @pytest.mark.skipif("not suitesparse") From 6a264344d708c4f0ee2ea155bad62f4545a19db5 Mon Sep 17 00:00:00 2001 From: Erik Welch Date: Thu, 12 Nov 2026 13:09:30 -0600 Subject: [PATCH 8/9] Compare in the type the dtypes unify to in Scalar.isclose without Numba Without Numba, Scalar.isclose passed the Python values straight to the isclose function. Every other path compares in the type the two dtypes unify to: Vector and Matrix isclose, and Scalar.isclose with Numba, which passes its values in the op's input types. INT64 and UINT64 unify to FP64, where 2**62 + 1 rounds to 2**62: Scalar.from_value(2**62 + 1).isclose( Scalar.from_value(2**62, dtype="UINT64"), rel_tol=0) with Numba, and Vector.isclose of the same: True Scalar.isclose without Numba: False, now True The values are converted to that type and back to Python numbers, so x - y still cannot wrap. --- graphblas/core/scalar.py | 9 ++++++--- graphblas/tests/test_scalar.py | 13 +++++++++++++ 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/graphblas/core/scalar.py b/graphblas/core/scalar.py index b3b78193d..0b711c258 100644 --- a/graphblas/core/scalar.py +++ b/graphblas/core/scalar.py @@ -320,8 +320,9 @@ def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False): # We can't yet call a UDF on a scalar as part of the spec, so let's do it ourselves isclose_func = isclose(rel_tol, abs_tol) if not _has_numba: - # Check if types are compatible - get_typed_op( + # Check if types are compatible, and compare in the type they unify to, as + # with Numba (INT64 with UINT64 is in FP64). Python ints do not wrap in x - y. + eq = get_typed_op( binary.eq, self.dtype, other.dtype, @@ -329,7 +330,9 @@ def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False): is_right_scalar=True, kind="binary", ) - return isclose_func(self.value, other.value) + x = eq.type.np_type.type(self.value).item() + y = eq.type2.np_type.type(other.value).item() + return isclose_func(x, y) isclose_func = get_typed_op( isclose_func, self.dtype, diff --git a/graphblas/tests/test_scalar.py b/graphblas/tests/test_scalar.py index d4168cd21..35439f34b 100644 --- a/graphblas/tests/test_scalar.py +++ b/graphblas/tests/test_scalar.py @@ -226,6 +226,19 @@ def test_isclose_integers_do_not_wrap(): assert left.isclose(right, rel_tol=rel_tol) is expected, (lt, lv, rt, rv) +def test_isclose_in_the_type_the_dtypes_unify_to(): + # INT64 and UINT64 unify to FP64, where 2**62 + 1 rounds to 2**62, so they are close + # even with no tolerance, as in Vector.isclose; with Numba and without + x = Scalar.from_value(2**62 + 1, dtype="INT64") + y = Scalar.from_value(2**62, dtype="UINT64") + assert x.isclose(y, rel_tol=0) + assert y.isclose(x, rel_tol=0) + v = Vector.from_coo([0], [2**62 + 1]) + assert v.isclose(Vector.from_coo([0], [2**62], "UINT64"), rel_tol=0) + # Two INT64 values are compared exactly + assert not x.isclose(Scalar.from_value(2**62, dtype="INT64"), rel_tol=0) + + def test_isclose_udt_raises(): """Reject ``isclose`` on UDTs with a clear message that points to isequal.""" arr_udt = dtypes.register_anonymous(np.dtype((np.int64, (3,))), "_IscloseUdt3") From 5047fa74f7d08ad73951028093c56f98fb3470f7 Mon Sep 17 00:00:00 2001 From: Erik Welch Date: Thu, 12 Nov 2026 13:09:33 -0600 Subject: [PATCH 9/9] Describe the coo indptr helpers by what they return, not by the old bug The indices_to_indptr docstring told how the old kernel lost entries, which the commit that fixed it already says. It now states the invariant: every empty row gets an empty range. The select test's docstring no longer names main as a reference point. --- graphblas/core/ss/matrix.py | 9 ++++----- graphblas/tests/test_matrix.py | 4 ++-- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/graphblas/core/ss/matrix.py b/graphblas/core/ss/matrix.py index 3205a1a69..d39eb9145 100644 --- a/graphblas/core/ss/matrix.py +++ b/graphblas/core/ss/matrix.py @@ -4564,11 +4564,10 @@ def issorted(arr): # pragma: no cover (numba) @njit def indices_to_indptr(indices, size): # pragma: no cover (numba) - """Calculate the indptr for e.g. CSR from sorted uint64 COO rows, each below size - 1. + """Calculate the indptr, e.g. for CSR, from uint64 COO rows, each below ``size - 1``. - Each row's count goes at ``indptr[row + 1]``, then a running sum. Setting only the - rows that appear left 0 for every empty row after the first one, so an empty row - in the middle moved entries to another row, and an empty last row lost them all. + Each row's count goes at ``indptr[row + 1]``, then a running sum makes the offsets, + so every empty row, wherever it is, gets an empty range. """ indptr = np.zeros(size, dtype=indices.dtype) one = np.uint64(1) @@ -4580,7 +4579,7 @@ def indices_to_indptr(indices, size): # pragma: no cover (numba) def _coo_indptr(indices, size, name, size_name): - """Return the indptr, of ``size + 1``, of sorted COO ``indices`` (rows or columns).""" + """Return the CSR (or CSC) indptr, of length ``size + 1``, for sorted COO rows (or columns).""" # Read 2-D indices as a flat buffer, as other imports do. The kernel would index # a 2-D array by its rows, past the end of its first axis. indices = ints_to_numpy_buffer(indices, np.uint64, name=name).ravel() diff --git a/graphblas/tests/test_matrix.py b/graphblas/tests/test_matrix.py index 2f10e5ed0..611e1a13e 100644 --- a/graphblas/tests/test_matrix.py +++ b/graphblas/tests/test_matrix.py @@ -1303,8 +1303,8 @@ def test_apply_indexunary(A): def test_select_by_a_transposed_comparison(): """``A.select(A.T == x)`` keeps the entries where ``A.T == x``, as a computed mask does. - It was rewritten to ``select.valueeq(A, x)``, which drops the transpose, and - for a matrix that is not square the DimensionMismatch (on main since #252). + Rewriting it to ``select.valueeq(A, x)`` would drop the transpose and, for a matrix + that is not square, the DimensionMismatch. """ S = Matrix.from_coo([0, 1], [1, 0], [1, 5], nrows=2, ncols=2) expected = Matrix.from_coo([1], [0], [5], nrows=2, ncols=2)