diff --git a/docs/user_guide/udt.rst b/docs/user_guide/udt.rst index cd2f93932..ea77dd82b 100644 --- a/docs/user_guide/udt.rst +++ b/docs/user_guide/udt.rst @@ -128,6 +128,8 @@ The following operators auto-lift to any UDT on first use: ``second``, ``any``, and ``pair`` work too, since they don't touch field values. - UnaryOps: ``ainv``, ``abs``. +- ``select`` with ``valueeq`` or ``valuene`` (``v.select("==", x)``), which + keeps the elements that ``==`` or ``!=`` finds. - Monoids: ``plus``, ``times``, ``min``, ``max``, ``any``. - Semirings combining a UDT-lifting monoid with a UDT-lifting BinaryOp (e.g., ``plus_times``, ``min_plus``, ``max_times``, ``any_first``). @@ -135,6 +137,11 @@ The following operators auto-lift to any UDT on first use: ``any_value``, ``count``. The positional ``agg.ss.first`` and ``agg.ss.last`` work on any dtype, including UDTs. +Except ``pair``, ``count``, ``agg.ss.first`` and ``agg.ss.last``, they need +Numba. Without it, the others raise ``KeyError`` on a UDT, and so do ``==`` and +``isequal`` on a Vector or Matrix of a UDT; a ``Scalar``'s ``==`` and ``isequal`` +compare in numpy. + For example: .. code-block:: python @@ -256,8 +263,9 @@ data. * - ...from an existing object: ``w << v``, ``v.dup(dtype=...)``, with a mask, an accumulator or ``replace`` - One pass through a cast op. - * - ...from an element-wise lifted op: ``apply``, ``ewise_mult``, - ``ewise_union``, and infix operators that use them + * - ...from an element-wise lifted op, or ``identity``, ``first`` or + ``second``: ``apply``, ``ewise_mult``, ``ewise_union``, and infix + operators that use them - The op writes the object's type itself, in the same call: ``int_udt_vec += 0.5`` truncates, as it does for an ``INT8`` vector. * - ...from anything else: ``mxm``, ``reduce`` into a Vector, ``ewise_add``, @@ -272,6 +280,13 @@ data. cast before it accumulates, as GraphBLAS converts a value of a built-in dtype, and so is a literal: an ``INT8[3]`` Scalar's ``s += 0.5`` adds 0, as an ``INT8`` Scalar's does. + * - A literal: ``v[0] = (1.5, 1)``, ``Scalar.from_value``, + ``Vector.from_scalar`` + - Converted as numpy assigns it to an element, so, unlike a literal that + an op makes an element, it need not fit: ``1.5`` in an integer field + truncates, text too long for its field is cut, and a datetime in a date + field loses its time of day. An int out of an integer field's range + raises ``OverflowError``, where an ``INT8`` Vector wraps it. * - Between a UDT and a built-in dtype, or layouts that do not correspond - ``DomainMismatch``. diff --git a/graphblas/core/dtypes.py b/graphblas/core/dtypes.py index 3ba87aff6..a5dddd017 100644 --- a/graphblas/core/dtypes.py +++ b/graphblas/core/dtypes.py @@ -112,10 +112,14 @@ def _get_udt_c_type_info(datatype): ``_gbudt_NNN`` name when the user didn't supply a name, or supplied one that isn't a valid C identifier or is too long for SuiteSparse (``udt_codegen._pick_c_type_name``), so the UDT can still take the JIT - path. Returns ``None`` only for the cases that remain unrepresentable in - C: an unsupported field type (object, datetime, string), or a field name - that collides with a C reserved word or isn't an identifier. Supports - record UDTs (including nested records) and array UDTs (e.g., ``FP64[3]``). + path. Returns ``None`` when the UDT has no C form here: a record laid out + other than as a C compiler lays out its struct (packed, as + ``np.dtype([...])`` is without ``align=True``), a field name that is not + a valid C identifier (a C reserved word, say), a field that is an array + (only an array UDT is flattened, to ``v[n]``), a field or element type + with no C type (object, datetime, text), or a UDT that is neither a + record nor an array. Supports record UDTs (including nested records) and + array UDTs (e.g., ``FP64[3]``). """ from .operator.udt_codegen import _udt_c_typedef @@ -212,6 +216,28 @@ def _numba_type(dtype): f"Numba, which python-graphblas uses for UDTs, cannot type an array of {element}, " f"as in {dtype}; use an array of {instead} instead" ) from exc + except numba.core.errors.NumbaError as exc: + # Numba's own error is "NumbaNotImplementedError: |V4", with terminal colors + raw = _raw_bytes_element(dtype) + if raw is None: # pragma: no cover (safety) + raise + raise ValueError( + f"Numba, which python-graphblas uses for UDTs, cannot type raw bytes ({raw}), " + f"as in {dtype}; use an array of uint8 instead, np.dtype((np.uint8, " + f"{raw.itemsize}))" + ) from exc + + +def _raw_bytes_element(dtype): + """The first void dtype with no fields (raw bytes) in dtype, or None.""" + if dtype.subdtype is not None: + return _raw_bytes_element(dtype.subdtype[0]) + if dtype.names is None: + return dtype if dtype.kind == "V" else None + for field_dtype, *_ in dtype.fields.values(): + if (raw := _raw_bytes_element(field_dtype)) is not None: + return raw + return None def register_new(name, dtype=None): @@ -605,13 +631,14 @@ def lookup_dtype(key, value=None): try: # Auto-register! return register_anonymous(key) - except Exception: - pass + except Exception as exc: + # Why it cannot be registered, such as an array of bool, which Numba cannot type. + reason = exc try: return lookup_dtype(key.literal_type) # For numba dtype inference except Exception: pass - raise ValueError(f"Unknown dtype: {key} of type {type(key)}") + raise ValueError(f"Unknown dtype: {key} of type {type(key)}") from reason def _raise_dtype_or_arraylike(cls_name, dtype, exc): diff --git a/graphblas/core/literals.py b/graphblas/core/literals.py index b9255093e..8945826ef 100644 --- a/graphblas/core/literals.py +++ b/graphblas/core/literals.py @@ -31,8 +31,8 @@ from .operator.monoid import Monoid from .udt import ( BUILTIN_UDT_BINARY_OPS, + _as_record_parts, _casts_safely_into, - _dict_to_record, _leaf_arrays, _promoted_record, _strip_leading_ones, @@ -102,7 +102,10 @@ def _literal_operand(dtype, value, op, *, left=False): return (math.inf if value > 0 else -math.inf), FP64 return value, literal_dtype if getattr(op, "_udt_types", True) is None: - # An op that does not take UDTs at all (``lt``) reports that itself. + # An op that does not take UDTs at all (``lt``) reports that itself, here, + # before the literal becomes an element, which can fail first (300 in an + # int8 field). + op[dtype] return value, dtype if isinstance(op, Monoid) and not op._anonymous: # Element-wise, a Monoid is its BinaryOp (``get_typed_op`` types it so). @@ -113,11 +116,14 @@ def _literal_operand(dtype, value, op, *, left=False): # truediv gives float leaves whatever the int, so, as a comparison, it # types an int out of an integer leaf's range exactly instead of raising. exact = is_comparison or op.name == "truediv" - if is_comparison and not isinstance(value, (np.generic, np.ndarray)): - np_type = dtype.np_type - parts = value - if isinstance(value, dict) and np_type.names is not None: - parts = _dict_to_record(np_type, value) + np_type = dtype.np_type + if is_comparison and np_type.kind == "S" and isinstance(value, str): + # A str, np.str_ included, that bytes cannot hold (numpy encodes ASCII) + # equals no element, as in numpy; one byte longer than any, it compares so. + if not value.isascii(): + value = b"\xff" * (np_type.itemsize + 1) + elif is_comparison and not isinstance(value, (np.generic, np.ndarray)): + parts = _as_record_parts(np_type, value) # First: making an int beyond its leaf infinity reshapes the part. _refuse_misfit_fields(dtype, np_type, parts) if _int_beyond_leaf(np_type, parts) is not None: @@ -311,13 +317,12 @@ def _literal_type(dtype, value, *, exact=False): """ value = _as_python_number(value) np_type = dtype.np_type - if isinstance(value, dict) and np_type.names is not None: - value = _dict_to_record(np_type, value) + value = _as_record_parts(np_type, value) _refuse_misfit_fields(dtype, np_type, value) if np_type.names is None and np_type.kind in "SU" and isinstance(value, (str, bytes)): if _keeps_value(value, np_type): return dtype - return lookup_dtype(np.dtype(f"{np_type.kind}{max(len(value), 1)}")) + return lookup_dtype(_own_text_type(np_type.kind, value)) if np_type.kind in "mM" and (date := _as_numpy_date(value, np_type)) is not None: if date.dtype.kind == np_type.kind and _keeps_value(date, np_type): return dtype @@ -339,6 +344,17 @@ def _literal_type(dtype, value, *, exact=False): return None +def _own_text_type(kind, text): + """Return the type of ``text`` that does not fit a text leaf of ``kind``: its own length. + + It keeps the leaf's kind, except that a str that bytes cannot hold (numpy + encodes ASCII) is a str, so beside bytes it equals nothing, as in numpy. + """ + if kind == "S" and isinstance(text, str) and not text.isascii(): + kind = "U" + return np.dtype(f"{kind}{max(len(text), 1)}") + + def _check_literal_fits(dtype, value): """Raise unless literal ``value`` may become an element of UDT ``dtype``, by its type. @@ -373,8 +389,7 @@ def _check_literal_fits(dtype, value): np_type = dtype.np_type if np_type.names is not None: target = np.zeros(1, dtype=np_type) - if isinstance(value, dict): - value = _dict_to_record(np_type, value) + value = _as_record_parts(np_type, value) else: # An array's elements, or the one value of a UDT that is neither (bytes). base, shape = np_type.subdtype or (np_type, ()) @@ -445,8 +460,7 @@ def _converts_exactly(dtype, value): np_type = dtype.np_type try: if np_type.names is not None: - if isinstance(value, dict): - value = _dict_to_record(np_type, value) + value = _as_record_parts(np_type, value) given = np.zeros(1, dtype=_object_layout(np_type)) converted = np.zeros(1, dtype=np_type) else: @@ -671,7 +685,7 @@ def leaf_dtype(leaf, part, *, in_record=False): # Text in a duration field, or text numpy cannot read as a date, # takes the field's type, and converting it reports what is wrong. if isinstance(part, (str, bytes)) and base.kind in "SU": - own = np.dtype(f"{base.kind}{max(len(part), 1)}") + own = _own_text_type(base.kind, part) elif (date := _as_numpy_date(part, base)) is not None: part = date own = date.dtype @@ -874,18 +888,19 @@ def _array_literal_udt(dtype, value): ``+ np.array([1, 2, 3])`` int64 ones. Its type is the array layout of its dtype: in the UDT's shape when it broadcasts to that (``np.array([1])`` beside ``INT8[3]``), which a literal that must become an element needs, - else in its own shape when the two broadcast together, so that an op's - result takes the shape numpy's would (``INT8[3]`` and a 2-by-3 array give - ``[2, 3]``). A 0-d structured array or ``np.void`` beside a record is its - record dtype. Anything else, such as a numeric array beside a record, is - converted into the UDT. + else in its own shape, so that an op's result takes the shape numpy's + would (``INT8[3]`` and a 2-by-3 array give ``[2, 3]``), and an op refuses + shapes that do not broadcast together, as for two UDTs. A 0-d structured + array or ``np.void`` beside a record is its record dtype. Anything else, + such as a numeric array beside a record, is converted into the UDT. """ np_type = dtype.np_type if value.dtype.names is not None: return lookup_dtype(value.dtype) if value.ndim == 0 else None # A bool array converts as it is (a BOOL array UDT cannot be registered), # and bool casts safely into every numeric element. - if np_type.names is not None or value.dtype.kind not in "iufc": + if np_type.names is not None or value.dtype.kind not in "iufc" or value.size == 0: + # An empty array's layout cannot be a UDT, so numpy says why it does not fit. return None shape = np_type.subdtype[1] if _strip_leading_ones(value.shape) == _strip_leading_ones(shape): @@ -894,9 +909,5 @@ def _array_literal_udt(dtype, value): try: np.broadcast_to(value, shape) except ValueError: - try: - np.broadcast_shapes(value.shape, shape) - except ValueError: - return None shape = value.shape return lookup_dtype(np.dtype((value.dtype, shape))) diff --git a/graphblas/core/matrix.py b/graphblas/core/matrix.py index 717a16ead..c14ba79ec 100644 --- a/graphblas/core/matrix.py +++ b/graphblas/core/matrix.py @@ -61,7 +61,15 @@ values_to_numpy_buffer, wrapdoc, ) -from .vector import Vector, VectorExpression, VectorIndexExpr, _isclose_recipe, _select_mask +from .vector import ( + Vector, + VectorExpression, + VectorIndexExpr, + _isclose_recipe, + _select_by_comparison, + _select_mask, + _udt_select_comparison, +) if backend == "suitesparse": from .ss.matrix import ss @@ -2809,9 +2817,10 @@ def select(self, op, thunk=None): if thunk is None: thunk = False # most basic form of 0 when unifying dtypes given = thunk + comparison = _udt_select_comparison(op, self.dtype) if type(thunk) is not Scalar: # A thunk beside a UDT becomes an element of it, so it must fit, as in apply. - thunk, dtype = _literal_operand(self.dtype, thunk, op) + thunk, dtype = _literal_operand(self.dtype, thunk, comparison or op) try: thunk = Scalar.from_value(thunk, dtype, is_cscalar=None, name="") except TypeError: @@ -2823,6 +2832,17 @@ def select(self, op, thunk=None): extra_message="Literal scalars also accepted.", op=op, ) + if comparison is not None: + mask_expr = self.apply(comparison, right=thunk) # typed now, so errors raise now + return MatrixExpression( + method_name, + None, + [self, thunk, _select_by_comparison, (self, mask_expr)], + expr_repr=f"{{0.name}}.select({op}, thunk={{1._expr_name}})", + nrows=self._nrows, + ncols=self._ncols, + dtype=self.dtype, + ) op = get_typed_op(op, self.dtype, thunk.dtype, is_right_scalar=True, kind="select") op, thunk = _value_select_in_own_type(op, self.dtype, thunk, given) self._expect_op(op, ("SelectOp", "IndexUnaryOp"), within=method_name, argname="op") diff --git a/graphblas/core/operator/base.py b/graphblas/core/operator/base.py index 0654c3d22..0192166a5 100644 --- a/graphblas/core/operator/base.py +++ b/graphblas/core/operator/base.py @@ -429,6 +429,10 @@ def __getitem__(self, type_): if type(type_) is tuple: from .utils import get_typed_op + if len(type_) != 2: + raise TypeError( + f"{self!r}[...] takes a dtype or a pair of dtypes, not a tuple of {len(type_)}" + ) dtype1, dtype2 = type_ dtype1 = lookup_dtype(dtype1) dtype2 = lookup_dtype(dtype2) diff --git a/graphblas/core/operator/binary.py b/graphblas/core/operator/binary.py index cb8abde73..dce4c7eda 100644 --- a/graphblas/core/operator/binary.py +++ b/graphblas/core/operator/binary.py @@ -846,7 +846,9 @@ def _compile_udt(self, dtype, dtype2): self.name, py_op, dtype, dtype2 ) elif self._numba_func is None: - raise KeyError(f"{self.name} does not work with {dtypes} types") + # Only without Numba does an op that lifts to UDTs get here. + needs = "; built-in ops work on UDTs only with Numba installed" if is_lifted else "" + raise KeyError(f"{self.name} does not work with {dtypes} types{needs}") else: numba_func = self._numba_func sig = (dtype.numba_type, dtype2.numba_type) diff --git a/graphblas/core/operator/indexbinary.py b/graphblas/core/operator/indexbinary.py index 41f7e3b3d..cf0f8f3a5 100644 --- a/graphblas/core/operator/indexbinary.py +++ b/graphblas/core/operator/indexbinary.py @@ -71,6 +71,17 @@ def _theta_to_scalar(theta, dtype): "(op[dtype](theta), or op(theta, dtype=...))." ) from exc raise + except OverflowError as exc: + # cffi's message, "int too big to convert", names neither the int nor the type. + if type(theta) is not int or dtype is not None and dtype._is_udt: + raise + if dtype is None: + raise OverflowError( + f"Python integer {theta} out of bounds for int64, the type of a Python int " + "theta; give theta as a Scalar of a type that holds it, or name the op's " + "type (op[dtype](theta), or op(theta, dtype=...))." + ) from exc + raise OverflowError(f"Python integer {theta} out of bounds for {dtype.np_type}") from exc def _refuse_empty_theta(theta): @@ -578,6 +589,11 @@ def _get_typed(self, dtype, dtype2, theta_type): def __getitem__(self, type_): if type(type_) is tuple and len(type_) == 3: return self._get_typed(*(lookup_dtype(dtype) for dtype in type_)) + if type(type_) is tuple and len(type_) != 2: + raise TypeError( + f"{self!r}[...] takes a dtype, or a tuple of two (x and y) or three (x, y and " + f"theta), not a tuple of {len(type_)}" + ) return super().__getitem__(type_) def _add(self, op, *, is_jit=False): diff --git a/graphblas/core/operator/monoid.py b/graphblas/core/operator/monoid.py index bae9c9db6..1880c0ed8 100644 --- a/graphblas/core/operator/monoid.py +++ b/graphblas/core/operator/monoid.py @@ -19,7 +19,7 @@ UINT64, lookup_dtype, ) -from ...exceptions import check_status_carg +from ...exceptions import EmptyObject, check_status_carg from .. import _has_numba, ffi, lib from ..utils import libget from .base import OpBase, ParameterizedUdf, TypedOpBase, _BinaryopJitDelegate, _hasop @@ -304,6 +304,12 @@ def _build(cls, name, binaryop, identity, *, is_idempotent=False, anonymous=Fals "a new BinaryOp (BinaryOp.register_new or BinaryOp.register_anonymous) " "and pass that." ) + from ..scalar import Scalar + + identities = identity.values() if isinstance(identity, Mapping) else [identity] + if any(isinstance(ident, Scalar) and ident._is_empty for ident in identities): + # A UDT monoid would get an identity of zero bytes, a wrong one for ``times``. + raise EmptyObject("identity is an empty Scalar; give a Scalar that holds a value") if name is None: name = binaryop.name new_type_obj = cls( @@ -378,7 +384,13 @@ def _compile_udt(self, dtype, dtype2): else: identity = Scalar.from_value(identity_val, dtype=ret_type, is_cscalar=True) new_monoid = ffi_new("GrB_Monoid*") - status = lib.GrB_Monoid_new_UDT(new_monoid, binaryop.gb_obj, identity.gb_obj) + if ret_type._is_udt: + status = lib.GrB_Monoid_new_UDT(new_monoid, binaryop.gb_obj, identity.gb_obj) + else: + # An is_udt op typed for a built-in type, as a monoid over a UDF built by + # ``_build`` is. GrB_Monoid_new_UDT refused it with an empty DomainMismatch. + func = libget(f"GrB_Monoid_new_{ret_type.name}") + status = func(new_monoid, binaryop.gb_obj, identity.gb_obj[0]) check_status_carg(status, "Monoid", new_monoid[0]) op = TypedUserMonoid( self, diff --git a/graphblas/core/operator/udf.py b/graphblas/core/operator/udf.py index 1fb644b73..850f5d6a8 100644 --- a/graphblas/core/operator/udf.py +++ b/graphblas/core/operator/udf.py @@ -58,17 +58,21 @@ def _numba_compile(numba_func, sig): also warn NumbaPendingDeprecationWarning that code may depend on this, and ``-W error`` (as in our tests) raises that warning instead of the error. The callers here handle both styles, so the warning only gets in the way. + Its text need not start with what the filter looks for: it may begin with + terminal colors or with "Failed in nopython mode pipeline", and a filter + matches from the start, hence the ``(?s).*``. """ # MAINT 2026-10-28: once Numba 0.61 is the minimum (0.60 made "new_style" the - # default; 0.61 removed "old_style"), drop this warning filter and the - # LoweringError branch in _compile_for_dtype. + # default; 0.61 removed "old_style"), drop this warning filter, the + # LoweringError branch in _compile_for_dtype, _OLD_STYLE_MISSING_FIELD_RE, and + # the unquoting in _summarize_numba_typing_error. if getattr(numba.core.config, "CAPTURED_ERRORS", None) != "old_style": numba_func.compile(sig) return with warnings.catch_warnings(): warnings.filterwarnings( "ignore", - "Code using Numba extension API maybe depending on 'old_style'", + r"(?s).*Code using Numba extension API maybe depending on 'old_style'", NumbaPendingDeprecationWarning, ) numba_func.compile(sig) @@ -144,18 +148,44 @@ def _compile_udf_for_udt(numba_func, sig, *, op_kind, op_name, dtypes): UnsupportedError, ...) so any compilation failure produces a single-line UDT diagnostic. Also catches the bare NotImplementedError Numba raises for an op it types but cannot lower, such as ``>`` (so - ``max``) on complex values. + ``max``) on complex values, and the bare KeyError for a missing field. """ + dtypes_str = ", ".join(str(d) for d in dtypes) try: _numba_compile(numba_func, sig) - except (NumbaError, NotImplementedError) as exc: - dtypes_str = ", ".join(str(d) for d in dtypes) - snippet = _summarize_numba_typing_error(exc) + except (NumbaError, NotImplementedError, KeyError) as exc: + if (field := _missing_field(exc)) is not None: + snippet = f"Field {field!r} was not found in the record" + else: + snippet = _summarize_numba_typing_error(exc) raise UdfParseError( f"{op_kind}.{op_name} does not work with ({dtypes_str}): {snippet}" ) from exc +# Before 0.60, Numba captures that KeyError "old_style": the TypingError it raises +# gives the KeyError (the quoted name) on the line after this one. +_OLD_STYLE_MISSING_FIELD_RE = re.compile( + r'Internal error at resolving type of attribute "(?P[^"]+)" of "[^"]*"\.' + r"\s*'(?P=name)'" +) + + +def _missing_field(exc): + """Return the field that ``x.name`` reads and the record lacks, if that is why ``exc``. + + Numba types ``x.name`` for a missing field with a bare KeyError of the name. + Before 0.60 a TypingError holds that KeyError instead. ``x["name"]`` gets + Numba's own "Field ... was not found" sentence: in a TypingError, or, in + 0.60, as a bare KeyError of it, which is not a name. + """ + if isinstance(exc, KeyError): + name = exc.args[0] if exc.args else None + return name if isinstance(name, str) and name.isidentifier() else None + match = _OLD_STYLE_MISSING_FIELD_RE.search(_ANSI_RE.sub("", str(exc))) + return None if match is None else match["name"] + + # Numba prefixes its actionable diagnostic line with one of these. _NUMBA_DIAG_PREFIXES = ( "No implementation of function", @@ -182,6 +212,10 @@ def _summarize_numba_typing_error(exc): """ lines = [_ANSI_RE.sub("", line).strip() for line in str(exc).splitlines()] for i, line in enumerate(lines): + # A KeyError, such as "Field ... was not found" for ``x["name"]``, prints as + # its quoted repr: bare in Numba 0.60, and inside the TypingError before it + if len(line) > 1 and line[0] == line[-1] == '"': + line = line[1:-1] if line.startswith(_NUMBA_DIAG_PREFIXES): if line.endswith(":"): for follow in lines[i + 1 :]: diff --git a/graphblas/core/operator/udt_codegen.py b/graphblas/core/operator/udt_codegen.py index 0ac54404a..99bd3e58a 100644 --- a/graphblas/core/operator/udt_codegen.py +++ b/graphblas/core/operator/udt_codegen.py @@ -33,6 +33,7 @@ _check_udt_pair, _get_udt_info, _iter_record_leaves, + _records_in_array_message, _udt_op_types, ) @@ -1120,6 +1121,10 @@ def compile_udt_unary_wrapper(op_name, py_op, dtype, out=None): f"Element-wise UDT ops require a record dtype (named fields) " f"or an array dtype (e.g., FP64[3])." ) + np_type = dtype.np_type + leaves = [d for *_, d in _iter_record_leaves(np_type)] if np_type.names else [np_type] + if any(leaf.base.names is not None for leaf in leaves): + raise KeyError(_records_in_array_message(f"unary.{op_name}", dtype)) kind, _detail = info if kind == "record": @@ -1160,20 +1165,27 @@ def fused_store_op(typed_op, out_type): and casts each element or leaf into ``out_type`` as it writes it, as a store casts (``udt._store_cast_op``). So a result goes into an object of another UDT type in the one GraphBLAS call that computes it, with no - temporary. Only the built-in ops that lift to UDTs have such a form; - ``udt._store_cast_op`` has already accepted the pair of types. + temporary. Only the built-in ops that lift to UDTs have such a form, + and ``identity``, ``first`` and ``second``, which store an operand as + it is; ``udt._store_cast_op`` has already accepted the pair of types. """ from .udf import _finalize_udt_op parent = getattr(typed_op, "parent", None) if parent is None or parent._anonymous: return None - if typed_op.opclass == "BinaryOp" and parent.name in BUILTIN_UDT_BINARY_OPS: + if typed_op.opclass == "UnaryOp" and parent.name == "identity": + from .unary import _udt_cast_op + + return _udt_cast_op(typed_op.type, out_type) + if typed_op.opclass == "BinaryOp" and ( + parent.name in BUILTIN_UDT_BINARY_OPS or parent.name in {"first", "second"} + ): from .binary import BinaryOp, TypedUserBinaryOp cls, typed_cls = BinaryOp, TypedUserBinaryOp dtype2 = typed_op.type2 - py_op = BUILTIN_UDT_BINARY_OPS[parent.name] + py_op = BUILTIN_UDT_BINARY_OPS.get(parent.name) elif typed_op.opclass == "UnaryOp" and parent.name in BUILTIN_UDT_UNARY_OPS: from .unary import TypedUserUnaryOp, UnaryOp @@ -1191,7 +1203,9 @@ def fused_store_op(typed_op, out_type): dtype = typed_op.type if (types := dtype if dtype2 is None else (dtype, dtype2)) in fused._udt_ops: return fused._udt_ops[types] - if dtype2 is None: + if py_op is None: + wrapper, sig = _pick_cast_wrapper(parent.name, dtype, dtype2, out_type) + elif dtype2 is None: wrapper, sig, _ = compile_udt_unary_wrapper(parent.name, py_op, dtype, out=out_type) else: wrapper, sig, _ = compile_udt_binary_wrapper( @@ -1199,6 +1213,26 @@ def fused_store_op(typed_op, out_type): ) return _finalize_udt_op(fused, dtype, dtype2, out_type, wrapper, sig, typed_cls) + def _pick_cast_wrapper(op_name, dtype, dtype2, out_type): + """Compile ``first`` or ``second`` storing the operand it picks as UDT ``out_type``. + + The operand is cast as a store casts it (:func:`compile_udt_cast_wrapper`); + the other is not read, whatever its type. Returns (wrapper_func, wrapper_sig). + """ + nt = numba.types + picked = "x_ptr" if op_name == "first" else "y_ptr" + cast, cast_sig = compile_udt_cast_wrapper(dtype if picked == "x_ptr" else dtype2, out_type) + wrapper = _compile_codegen( + f"def _pick(z_ptr, x_ptr, y_ptr):\n cast(z_ptr, {picked})\n", + func_name="_pick", + source_label=f"", + extra_ns={"cast": numba.njit(cast)}, + ) + z_type, picked_type = cast_sig.args + if picked == "x_ptr": + return wrapper, nt.void(z_type, picked_type, nt.voidptr) + return wrapper, nt.void(z_type, nt.voidptr, picked_type) + # JIT C code generators below. @@ -1450,13 +1484,15 @@ def _c_expr_unary(py_op, operand, field_dtype): raise ValueError(f"Unknown unary C op: {py_op}") # pragma: no cover (safety) -def _make_jit_c_definition(op_name, py_op, dtype, arity): +def _make_jit_c_definition(op_name, py_op, dtype, arity, out=None): """Generate a JIT C function definition for a UDT operator. Returns ``(c_name, c_defn)``, or ``None`` if the dtype can't be expressed in C (unsupported field types, or names that collide with C reserved words). An op not defined on the field types (``min`` on a complex field) never - gets here: ``udt._element_types`` rejects it when the op is typed. + gets here: ``udt._element_types`` rejects it when the op is typed. ``out`` + is the result type when it is not ``dtype`` (``abs`` of complex leaves), + which SuiteSparse must already know by its C name. """ np_type = dtype.np_type # Prefer the C name SuiteSparse already has for this type. ``GxB_JIT_C_NAME`` @@ -1471,11 +1507,14 @@ def _make_jit_c_definition(op_name, py_op, dtype, arity): return None type_name, _typedef = typedef_info c_name = f"{op_name}_{type_name}" + z_type_name = type_name if out is None else out.jit_c_name + if z_type_name is None: + return None if arity == 2: - params = f"{type_name} *z, const {type_name} *x, const {type_name} *y" + params = f"{z_type_name} *z, const {type_name} *x, const {type_name} *y" else: - params = f"{type_name} *z, const {type_name} *x" + params = f"{z_type_name} *z, const {type_name} *x" if np_type.names is not None: # Use leaf C paths so nested records emit ``z->outer.inner_a = ...`` @@ -1740,7 +1779,7 @@ def set_jit_c_mixed_on_op(gb_obj, op_name, py_op, dtype, dtype2, ret_type, input return result -def set_jit_c_on_op(gb_obj, op_name, py_op, dtype, set_string_func, arity=2): +def set_jit_c_on_op(gb_obj, op_name, py_op, dtype, set_string_func, arity=2, out=None): """Generate and set JIT C name+definition on a GrB operator, if possible. Returns ``(c_name, c_defn)`` when the dtype is expressible in C, ``None`` @@ -1749,7 +1788,7 @@ def set_jit_c_on_op(gb_obj, op_name, py_op, dtype, set_string_func, arity=2): Callers check ``_has_jit_set`` first: without it, SuiteSparse has no ``set_string_func`` to pass. """ - result = _make_jit_c_definition(op_name, py_op, dtype, arity) + result = _make_jit_c_definition(op_name, py_op, dtype, arity, out) if result is None: return None if not _set_jit_c_strings(gb_obj, *result, set_string_func): diff --git a/graphblas/core/operator/unary.py b/graphblas/core/operator/unary.py index c879abe80..53ab96343 100644 --- a/graphblas/core/operator/unary.py +++ b/graphblas/core/operator/unary.py @@ -123,6 +123,7 @@ def _one(x): if _has_numba: from ..udt import BUILTIN_UDT_UNARY_OPS as _BUILTIN_UDT_UNARY_OPS + from ..udt import _unary_op_type from .udt_codegen import ( _has_jit_set, _maybe_warn_jit_skipped, @@ -287,8 +288,11 @@ def _compile_udt(self, dtype, dtype2): # Auto-generate an element-wise UDT function for the built-in # unary ops (``ainv``, ``abs``). py_op = _BUILTIN_UDT_UNARY_OPS[self.name] + # ``abs`` on a complex leaf gives a real one, written as it computes. + out = _unary_op_type(self.name, dtype) + out = None if out is dtype else out unary_wrapper, wrapper_sig, ret_type = compile_udt_unary_wrapper( - self.name, py_op, dtype + self.name, py_op, dtype, out=out ) op = _finalize_udt_op( self, dtype, None, ret_type, unary_wrapper, wrapper_sig, TypedUserUnaryOp @@ -298,7 +302,7 @@ def _compile_udt(self, dtype, dtype2): # there, so it must not be evaluated as an argument first. if _has_jit_set: op._jit_c_info = set_jit_c_on_op( - op.gb_obj, self.name, py_op, dtype, lib.GrB_UnaryOp_set_String, arity=1 + op.gb_obj, self.name, py_op, dtype, lib.GrB_UnaryOp_set_String, 1, out ) _maybe_warn_jit_skipped(op._jit_c_info, self.name, dtype.name) return op diff --git a/graphblas/core/operator/utils.py b/graphblas/core/operator/utils.py index e8413a168..1433d2d4d 100644 --- a/graphblas/core/operator/utils.py +++ b/graphblas/core/operator/utils.py @@ -18,6 +18,7 @@ from ...dtypes import UINT64, lookup_dtype, unify from ..expr import InfixExprBase from ..udt import _monoid_pair_type +from ..utils import _Pointer from .base import ( _SS_OPERATORS, OpBase, @@ -364,6 +365,8 @@ def _resolve_index_expr(ns, modname, expr, callname, opname): _check_index_expr(modname, expr, callname, opname) tensor = expr.args[0] thunk = expr.args[1] + if isinstance(thunk, _Pointer): # A UDT Scalar + thunk = thunk.val if expr._is_scalar: # Handle ScalarExpressions that change their arguments to Vector, before the # bound below reads the thunk's value. A comparison (an ewise_mult) has both diff --git a/graphblas/core/scalar.py b/graphblas/core/scalar.py index af3681e4e..cabbde48c 100644 --- a/graphblas/core/scalar.py +++ b/graphblas/core/scalar.py @@ -242,6 +242,9 @@ def isequal(self, other, *, check_dtype=False): if _literal_shape_differs(self.dtype, other): return False dtype = _literal_type(self.dtype, other, exact=True) or self.dtype + elif type(other) is int and not -(2**63) <= other < 2**63: + # INT64, a Python int's dtype, cannot hold it; compare exactly, as ``==`` does. + return not self._is_empty and self.value == other try: other = Scalar.from_value(other, dtype, is_cscalar=None, name="s_isequal") except TypeError: @@ -299,6 +302,14 @@ def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False): if type(other) is not Scalar: if other is None: return self._is_empty + if type(other) is int and not -(2**63) <= other < 2**63: + # INT64, a Python int's dtype, cannot hold it; compare in Python, as isequal does + if self._is_empty: + return False + try: + return _isclose(rel_tol, abs_tol)(self.value, other) + except OverflowError: # An int too large for a float, beside a float + return False try: other = Scalar.from_value(other, is_cscalar=None, name="s_isclose") except TypeError: diff --git a/graphblas/core/ss/_utils.py b/graphblas/core/ss/_utils.py index fafa06107..e68f3d589 100644 --- a/graphblas/core/ss/_utils.py +++ b/graphblas/core/ss/_utils.py @@ -4,9 +4,10 @@ from suitesparse_graphblas.utils import unclaim_buffer +from ... import dtypes from ...exceptions import _error_code_lookup from .. import ffi, lib -from ..dtypes import _string_to_dtype +from ..dtypes import DataType, _string_to_dtype def _unclaim_taken(*arrays_and_pointers): @@ -52,7 +53,8 @@ def _resolve_serialized_dtype(data_obj, data_nbytes, obj_kind): short name may not match any Python-side registration (e.g., when an anonymous UDT was registered by-shape on the producer side). Fall back to ``GrB_NAME`` (the numpy repr that ``serialize`` writes into the - object) when the short name fails to resolve. + object) when the short name fails to resolve, and then to a type + registered from C source under the short name. ``obj_kind`` is ``"Matrix"`` or ``"Vector"`` and only shows up in error messages. @@ -68,7 +70,8 @@ def _resolve_serialized_dtype(data_obj, data_nbytes, obj_kind): return _string_to_dtype(dtype_name) except (ValueError, TypeError, SyntaxError) as exc: orig_error = exc - if hasattr(lib, "GxB_Serialized_get_String"): + has_grb_name = hasattr(lib, "GxB_Serialized_get_String") + if has_grb_name: dtype_size = ffi.new("size_t*") info = lib.GxB_Serialized_get_SIZE(data_obj, dtype_size, lib.GrB_NAME, data_nbytes) if info != lib.GrB_SUCCESS: @@ -78,12 +81,21 @@ def _resolve_serialized_dtype(data_obj, data_nbytes, obj_kind): # This checks the blob as GxB_Serialized_get_SIZE just did, so it cannot fail now. if info != lib.GrB_SUCCESS: # pragma: no cover (safety) raise _error_code_lookup[info](f"{obj_kind} deserialize failed to get the name") - return _string_to_dtype(ffi.string(dtype_char).decode()) - # No GrB_NAME fallback available (older SS); surface the original - # resolution failure with the actual dtype name attached. + try: + return _string_to_dtype(ffi.string(dtype_char).decode()) + except (ValueError, TypeError, SyntaxError) as exc: + orig_error = exc + # A type registered from C source (``dtypes.ss.register_new``) has its own name as its + # C name, and no numpy repr in GrB_NAME (literal_eval("") raised a bare SyntaxError). + # GraphBLAS checks that it fits the blob. + c_source_type = getattr(dtypes.ss, dtype_name, None) if dtype_name else None + if type(c_source_type) is DataType: + return c_source_type + # Such as a type from C source that this process has not registered, or a UDT whose + # numpy repr GrB_NAME could not keep (or, before SuiteSparse 9, has no GrB_NAME) raise ValueError( - f"{obj_kind} deserialize cannot resolve dtype {dtype_name!r} and " - f"GxB_Serialized_get_String is unavailable in this SuiteSparse build" + f"{obj_kind} deserialize cannot tell the dtype from the type name " + f"{dtype_name!r}; pass it as dtype=" ) from orig_error diff --git a/graphblas/core/ss/binary.py b/graphblas/core/ss/binary.py index d940ac342..8cd9f303b 100644 --- a/graphblas/core/ss/binary.py +++ b/graphblas/core/ss/binary.py @@ -2,35 +2,12 @@ from ...dtypes import lookup_dtype from ...exceptions import check_status_carg from .. import NULL, ffi, lib -from ..operator.base import TypedOpBase -from ..operator.binary import BinaryOp, TypedUserBinaryOp from . import _IS_SSGB7 from ._utils import _jit_c_source ffi_new = ffi.new -class TypedJitBinaryOp(TypedOpBase): - __slots__ = "_monoid", "_jit_c_definition" - opclass = "BinaryOp" - - def __init__(self, parent, name, type_, return_type, gb_obj, jit_c_definition, dtype2=None): - super().__init__(parent, name, type_, return_type, gb_obj, name, dtype2=dtype2) - self._monoid = None - self._jit_c_definition = jit_c_definition - - @property - def jit_c_definition(self): - return self._jit_c_definition - - monoid = TypedUserBinaryOp.monoid - commutes_to = TypedUserBinaryOp.commutes_to - _semiring_commutes_to = TypedUserBinaryOp._semiring_commutes_to - is_commutative = TypedUserBinaryOp.is_commutative - type2 = TypedUserBinaryOp.type2 - __call__ = TypedUserBinaryOp.__call__ - - def register_new(name, jit_c_definition, left_type, right_type, ret_type): """Register a new BinaryOp using the SuiteSparse:GraphBLAS JIT compiler. @@ -139,3 +116,31 @@ def register_new(name, jit_c_definition, left_type, right_type, ret_type): rv._add(op, is_jit=True) setattr(module, funcname, rv) return rv + + +# Imported below register_new: importing the operator package initializes the +# operators, which imports graphblas.binary.ss, which takes register_new from here. +# Imported at the top, they made this module fail as a process's first import. +from ..operator.base import TypedOpBase # noqa: E402 isort:skip +from ..operator.binary import BinaryOp, TypedUserBinaryOp # noqa: E402 isort:skip + + +class TypedJitBinaryOp(TypedOpBase): + __slots__ = "_monoid", "_jit_c_definition" + opclass = "BinaryOp" + + def __init__(self, parent, name, type_, return_type, gb_obj, jit_c_definition, dtype2=None): + super().__init__(parent, name, type_, return_type, gb_obj, name, dtype2=dtype2) + self._monoid = None + self._jit_c_definition = jit_c_definition + + @property + def jit_c_definition(self): + return self._jit_c_definition + + monoid = TypedUserBinaryOp.monoid + commutes_to = TypedUserBinaryOp.commutes_to + _semiring_commutes_to = TypedUserBinaryOp._semiring_commutes_to + is_commutative = TypedUserBinaryOp.is_commutative + type2 = TypedUserBinaryOp.type2 + __call__ = TypedUserBinaryOp.__call__ diff --git a/graphblas/core/ss/matrix.py b/graphblas/core/ss/matrix.py index 9ad43721a..19d9963a7 100644 --- a/graphblas/core/ss/matrix.py +++ b/graphblas/core/ss/matrix.py @@ -23,6 +23,7 @@ _CArray, _distinct_buffers, _MatrixArray, + _ownable_copy, _Pointer, check_values_size, get_order, @@ -63,6 +64,69 @@ def get_udt(it): return get_udt +def _claim_values_2d(cdata, cdata_size, nrows, ncols, dtype, is_c_order): + """Like ``claim_buffer_2d``, for values, but keep each array-UDT element contiguous. + + GraphBLAS stores an array UDT's element as one contiguous block. Given a subarray + dtype in column order, ``claim_buffer_2d`` made the element's own axes the + slowest, so ``values[i, j]`` was not the entry at ``(i, j)``. Given more than + ``nrows * ncols`` elements, it raised ValueError for any subarray dtype. GraphBLAS + may hand over more: SuiteSparse 7's memory pool rounds an allocation up to a + power of two, and a full matrix resized smaller keeps its buffer. Row order with + an exactly sized buffer still goes to it: its array owns the buffer, so packing + the values back can take them without a copy. + """ + if not dtype.shape or (is_c_order and cdata_size == nrows * ncols): + return claim_buffer_2d(ffi, cdata, cdata_size, nrows, ncols, dtype, is_c_order) + values = claim_buffer(ffi, cdata, cdata_size, dtype)[: nrows * ncols] + if is_c_order: + return values.reshape(nrows, ncols, *dtype.shape) + return values.reshape(ncols, nrows, *dtype.shape).swapaxes(0, 1) + + +def _values_in_column_order(values, dtype, copy): + """Return ``(values, dtype, owner, c_contiguous)`` for a fullc or bitmapc import. + + GraphBLAS takes ``owner``'s memory, which holds the elements in column order: as + ``owner`` if ``c_contiguous``, else as ``owner.T``. An array UDT keeps each + element contiguous. Given one element per entry (``values[i, j]``), the two entry + axes are swapped into a copy; given a list of elements already in column order, + as ``raw=True`` exports them, numbers flat in that order, or one element (for an + iso value), they are taken as they are. Fortran order of the whole array, as for + other types, put the element's own axes slowest, and GraphBLAS read each element + from scattered numbers. + """ + # Converted once, to see the shape (and infer an array UDT from it). An array UDT's + # layout is decided below; other values go straight to Fortran order. + given = values + if isinstance(given, np.ndarray): + array_udt = given.ndim > 2 if dtype is None else bool(lookup_dtype(dtype).np_type.shape) + else: + array_udt = False + values, dtype = values_to_numpy_buffer( + values, dtype, order="K" if array_udt else "F", subarray_after=2 + ) + # A conversion that made a new array (from a list, or to another dtype or order) + # needs no second copy to be safe to hand over + if copy and not (isinstance(given, np.ndarray) and np.may_share_memory(values, given)): + copy = False + sub = dtype.np_type.shape + if sub and values.ndim == 2 + len(sub) and values.shape[2:] == sub: + return values, dtype, _ownable_copy(values.swapaxes(0, 1), "C"), True + if sub and values.ndim > 1 and values.shape[1:] != sub and values.shape != sub: + dims = ", ".join(map(str, sub)) + raise ValueError( + f"values of {dtype} in column order must have shape (nrows, ncols, {dims}), " + f"an element per entry, or (n, {dims}) or (n,), listed in column order, or " + f"({dims}), one element; " + f"got {values.shape}" + ) + values, dtype = values_to_numpy_buffer( + values, dtype, copy=copy, ownable=True, order="C" if sub else "F", subarray_after=2 + ) + return values, dtype, values, bool(sub) + + def head(matrix, n=10, dtype=None, *, sort=False): """Like ``matrix.to_coo()``, but only returns the first n elements. @@ -620,6 +684,9 @@ def export(self, format=None, *, sort=False, give_ownership=False, raw=False, ** Return values - Note: for ``raw=True``, arrays may be larger than specified. + - Note: for an array UDT, ``values`` has the element's own axes after those + shown, so with ``raw=False`` ``values[i, j]`` is the element at ``(i, j)``. + In "bitmapc" and "fullc" it is then not Fortran-ordered. - "csr" format - indptr : ndarray(dtype=uint64, ndim=1, size=nrows + 1) - col_indices : ndarray(dtype=uint64, ndim=1, size=nvals) @@ -814,18 +881,23 @@ def _export(self, format=None, *, sort=False, give_ownership=False, raw=False, m columns = _CArray(size=nvals, name="&columns_array") scalar = _scalar_index("s_nvals") scalar.value = nvals + dtype_name = "UDT" if parent.dtype._is_udt else parent.dtype.name call( - f"GrB_Matrix_extractTuples_{parent.dtype.name}", + f"GrB_Matrix_extractTuples_{dtype_name}", [rows, columns, None, _Pointer(scalar), parent], ) - value = parent.reduce_scalar(monoid.any, allow_empty=False).new().value + # One value, of shape (1, *element) for an array UDT. It is read with + # an iterator, which needs no Numba for a UDT, as reduce did. + values = np.zeros(1, dtype) # an iso Matrix with no entries has no value + if nvals: + values[0] = self.iso_value.value rv = { "format": "coo", "nrows": nrows, "ncols": ncols, "rows": rows.array, "cols": columns.array, - "values": np.array([value], dtype=dtype), + "values": values, "is_iso": True, } else: @@ -1140,14 +1212,8 @@ def _export(self, format=None, *, sort=False, give_ownership=False, raw=False, m if values.size > 1: # pragma: no cover (suitesparse) values = values[:1] else: - values = claim_buffer_2d( - ffi, - Ax[0], - Ax_size[0] // dtype.itemsize, - nrows, - ncols, - dtype, - is_c_order, + values = _claim_values_2d( + Ax[0], Ax_size[0] // dtype.itemsize, nrows, ncols, dtype, is_c_order ) rv = {"bitmap": bitmap, "nvals": nvals_[0]} if raw: @@ -1186,8 +1252,8 @@ def _export(self, format=None, *, sort=False, give_ownership=False, raw=False, m rv = {"nrows": nrows, "ncols": ncols} else: is_c_order = format == "fullr" - values = claim_buffer_2d( - ffi, Ax[0], Ax_size[0] // dtype.itemsize, nrows, ncols, dtype, is_c_order + values = _claim_values_2d( + Ax[0], Ax_size[0] // dtype.itemsize, nrows, ncols, dtype, is_c_order ) rv = {} else: @@ -2343,7 +2409,9 @@ def import_bitmapc( If True, perform a zero-copy data transfer from input numpy arrays to GraphBLAS if possible. To give ownership of the underlying memory buffers to GraphBLAS, the arrays must: - - be FORTRAN contiguous + - be FORTRAN contiguous, or for an array UDT, be a C-contiguous list + of elements in column order, as ``raw=True`` exports them (values + of shape ``(nrows, ncols, *element)`` are always copied) - have the correct dtype (bool for bitmap) - own its own data, allocated by GraphBLAS (arrays from ``export``, ``unpack`` and ``suitesparse_graphblas.utils.empty``) or by NumPy's @@ -2462,16 +2530,14 @@ def _import_bitmapc( ) if method == "pack": dtype = matrix.dtype - values, dtype = values_to_numpy_buffer( - values, dtype, copy=copy, ownable=True, order="F", subarray_after=2 - ) - bitmap, values = _distinct_buffers(bitmap, values) + values, dtype, owner, c_contiguous = _values_in_column_order(values, dtype, copy) + bitmap, owner = _distinct_buffers(bitmap, owner) if method == "import": nrows, ncols = get_shape(nrows, ncols, dtype, bitmap=bitmap, values=values) else: nrows, ncols = matrix.shape Ab = ffi_new("int8_t**", ffi.from_buffer("int8_t*", bitmap.T)) - Ax = ffi_new("void**", ffi.from_buffer("void*", values.T)) + Ax = ffi_new("void**", ffi.from_buffer("void*", owner if c_contiguous else owner.T)) if nvals is None: if bitmap.size == nrows * ncols: nvals = np.count_nonzero(bitmap) @@ -2491,12 +2557,12 @@ def _import_bitmapc( Ab, Ax, bitmap.nbytes, - values.nbytes, + owner.nbytes, is_iso, nvals, NULL if desc is None else desc._carg, ) - _unclaim_taken((bitmap, Ab), (values, Ax)) + _unclaim_taken((bitmap, Ab), (owner, Ax)) if method == "import": check_status_carg( status, @@ -2718,7 +2784,9 @@ def import_fullc( If True, perform a zero-copy data transfer from input numpy arrays to GraphBLAS if possible. To give ownership of the underlying memory buffers to GraphBLAS, the arrays must: - - be FORTRAN contiguous + - be FORTRAN contiguous, or for an array UDT, be a C-contiguous list + of elements in column order, as ``raw=True`` exports them (values + of shape ``(nrows, ncols, *element)`` are always copied) - have the correct dtype - own its own data, allocated by GraphBLAS (arrays from ``export``, ``unpack`` and ``suitesparse_graphblas.utils.empty``) or by NumPy's @@ -2824,14 +2892,12 @@ def _import_fullc( copy = not take_ownership or secure_import if method == "pack": dtype = matrix.dtype - values, dtype = values_to_numpy_buffer( - values, dtype, copy=copy, order="F", ownable=True, subarray_after=2 - ) + values, dtype, owner, c_contiguous = _values_in_column_order(values, dtype, copy) if method == "import": nrows, ncols = get_shape(nrows, ncols, dtype, values=values) else: nrows, ncols = matrix.shape - Ax = ffi_new("void**", ffi.from_buffer("void*", values.T)) + Ax = ffi_new("void**", ffi.from_buffer("void*", owner if c_contiguous else owner.T)) if method == "import": mhandle = ffi_new("GrB_Matrix*") args = (dtype._carg, nrows, ncols) @@ -2841,9 +2907,9 @@ def _import_fullc( check_values_size(values, dtype, nrows * ncols, is_iso) desc = get_descriptor(secure_import=secure_import, **opts) status = getattr(lib, f"GxB_Matrix_{method}_FullC")( - mhandle, *args, Ax, values.nbytes, is_iso, NULL if desc is None else desc._carg + mhandle, *args, Ax, owner.nbytes, is_iso, NULL if desc is None else desc._carg ) - _unclaim_taken((values, Ax)) + _unclaim_taken((owner, Ax)) if method == "import": check_status_carg( status, diff --git a/graphblas/core/ss/unary.py b/graphblas/core/ss/unary.py index 695dec5c1..442742aca 100644 --- a/graphblas/core/ss/unary.py +++ b/graphblas/core/ss/unary.py @@ -2,29 +2,12 @@ from ...dtypes import lookup_dtype from ...exceptions import check_status_carg from .. import NULL, ffi, lib -from ..operator.base import TypedOpBase -from ..operator.unary import TypedUserUnaryOp, UnaryOp from . import _IS_SSGB7 from ._utils import _jit_c_source ffi_new = ffi.new -class TypedJitUnaryOp(TypedOpBase): - __slots__ = "_jit_c_definition" - opclass = "UnaryOp" - - def __init__(self, parent, name, type_, return_type, gb_obj, jit_c_definition): - super().__init__(parent, name, type_, return_type, gb_obj, name) - self._jit_c_definition = jit_c_definition - - @property - def jit_c_definition(self): - return self._jit_c_definition - - __call__ = TypedUserUnaryOp.__call__ - - def register_new(name, jit_c_definition, input_type, ret_type): """Register a new UnaryOp using the SuiteSparse:GraphBLAS JIT compiler. @@ -120,3 +103,25 @@ def register_new(name, jit_c_definition, input_type, ret_type): rv._add(op, is_jit=True) setattr(module, funcname, rv) return rv + + +# Imported below register_new: importing the operator package initializes the +# operators, which imports graphblas.unary.ss, which takes register_new from here. +# Imported at the top, they made this module fail as a process's first import. +from ..operator.base import TypedOpBase # noqa: E402 isort:skip +from ..operator.unary import TypedUserUnaryOp, UnaryOp # noqa: E402 isort:skip + + +class TypedJitUnaryOp(TypedOpBase): + __slots__ = "_jit_c_definition" + opclass = "UnaryOp" + + def __init__(self, parent, name, type_, return_type, gb_obj, jit_c_definition): + super().__init__(parent, name, type_, return_type, gb_obj, name) + self._jit_c_definition = jit_c_definition + + @property + def jit_c_definition(self): + return self._jit_c_definition + + __call__ = TypedUserUnaryOp.__call__ diff --git a/graphblas/core/udt.py b/graphblas/core/udt.py index 27ebc7b90..0f91d7966 100644 --- a/graphblas/core/udt.py +++ b/graphblas/core/udt.py @@ -218,6 +218,13 @@ def _complex_op_error(op_name, dtype, dtype2): ) +def _records_in_array_message(op_label, dtypes_str): + return ( + f"{op_label} does not work with ({dtypes_str}): a built-in op lifts to the numbers " + "in an array or a record, not to records in an array; register a custom op for it." + ) + + def _element_types(op_name, dtype, dtype2, elem_x, elem_y): """Return ``(in_x, in_y, out)``, the numpy dtypes ``binary.{op_name}`` uses on one element. @@ -237,6 +244,9 @@ def _element_types(op_name, dtype, dtype2, elem_x, elem_y): # down with ``NotImplementedError: No definition for lowering lt``. if op_name in _OPS_NOT_FOR_COMPLEX and "c" in (elem_x.kind, elem_y.kind): raise _complex_op_error(op_name, dtype, dtype2) + if elem_x.names is not None or elem_y.names is not None: + # Records in an array (an array UDT or field); Numba cannot add two records. + raise KeyError(_records_in_array_message(f"binary.{op_name}", f"{dtype}, {dtype2}")) builtin_x = _builtin_dtype(elem_x) builtin_y = _builtin_dtype(elem_y) if builtin_x is None or builtin_y is None: @@ -315,6 +325,33 @@ def _monoid_pair_type(binaryop_name, dtype, dtype2): return _udt_op_types(op_name, dtype, dtype2)[0] +def _unary_op_type(op_name, dtype): + """Return the type ``unary.{op_name}`` gives on UDT ``dtype``: each leaf as on its built-in. + + Of the ops that lift, only ``abs`` changes a leaf's type: a complex number's + magnitude is real, as ``unary.abs`` on FC64 gives FP64 (and numpy's ``abs`` + a float), also where the backend has no complex types. + """ + if op_name != "abs": + return dtype + + def leaf_type(leaf): + base, shape = leaf.subdtype or (leaf, ()) + if base.kind == "c": + base = np.empty(0, base).real.dtype + return np.dtype((base, shape)) if shape else base + + np_type = dtype.np_type + if np_type.names is None: + result = leaf_type(np_type) + return dtype if result == np_type else lookup_dtype(result) + leaves = [leaf for _py, _c, leaf in _iter_record_leaves(np_type)] + result_leaves = [leaf_type(leaf) for leaf in leaves] + if result_leaves == leaves: + return dtype + return lookup_dtype(_promoted_record([np_type], iter(result_leaves))) + + def _promoted_record(records, result_leaves): """Build a record laid out like ``records`` with ``result_leaves`` for its leaves. @@ -489,8 +526,9 @@ def _fused_store_expr(expr, output_dtype): casts each element as it writes it (``udt_codegen.fused_store_op``), so no temporary of the result's own type is made: data GraphBLAS works on is often too large to hold twice. That works for an element-wise - computation by a built-in op that lifts to UDTs (``apply``, - ``ewise_mult``, ``ewise_union``). ``ewise_add`` copies an entry present in + computation by a built-in op that lifts to UDTs, or by ``identity``, + ``first`` or ``second`` (``apply``, ``ewise_mult``, ``ewise_union``); + ``any`` may pick either operand. ``ewise_add`` copies an entry present in only one input without the op, and a semiring or reduction would add in the narrower type, so those and any other computation raise. """ @@ -509,10 +547,10 @@ def _fused_store_expr(expr, output_dtype): raise DomainMismatch( f"cannot store a result of type {expr.dtype} in an object of type " f"{output_dtype}: GraphBLAS cannot cast a UDT, and only a built-in op that lifts " - "to UDTs, in apply, ewise_mult or ewise_union, can write another type as it " - f"computes ({what} cannot), so this would need a temporary copy of the whole " - "result. Make it explicitly by storing the computed result, as in " - "`obj << expr.new()`, which casts it in one more pass." + "to UDTs, or identity, first or second, in apply, ewise_mult or ewise_union, can " + f"write another type as it computes ({what} cannot), so this would need a " + "temporary copy of the whole result. Make it explicitly by storing the computed " + "result, as in `obj << expr.new()`, which casts it in one more pass." ) rv = copy.copy(expr) rv.op, rv.cfunc_name, rv.dtype = fused, cfunc_name, output_dtype @@ -785,8 +823,7 @@ def _element_bytes(dtype, value): # as junk however it writes the record. if np_type.subdtype is None: arr = np.zeros(1, dtype=np_type) - if isinstance(value, dict) and np_type.names is not None: - value = _dict_to_record(np_type, value) + value = _as_record_parts(np_type, value) else: arr = np.zeros(np_type.subdtype[1], dtype=np_type.subdtype[0]) value = _nats_in_leaf_units(_view_if_same_layout(value, arr.dtype), arr.dtype) @@ -818,16 +855,27 @@ def _element_value(dtype, data): return rv.view(base).reshape(shape) -def _dict_to_record(np_type, d): - """Converts e.g. ``{"x": 1, "y": 2.3}`` to ``(1, 2.3)``.""" - rv = [] - for name, (dtype, _) in np_type.fields.items(): - val = d[name] - if dtype.names is not None and isinstance(val, dict): - rv.append(_dict_to_record(dtype, val)) - else: - rv.append(val) - return tuple(rv) +def _as_record_parts(np_type, value): + """Return ``value`` for an element of ``np_type`` with each record in it as a tuple. + + A record literal may be a tuple, list or dict laid out like the element, but + numpy takes only a tuple: ``{"x": 1, "y": 2.3}`` and ``[1, 2.3]`` become + ``(1, 2.3)``, at every level of a nested record. A list of another length + becomes a tuple too, which numpy refuses by its length; an array field + keeps its list. + """ + if np_type.names is None: + return value + if isinstance(value, dict): + value = [value[name] for name in np_type.names] + elif isinstance(value, list) and len(value) != len(np_type.names): + return tuple(value) + elif not isinstance(value, (list, tuple)) or len(value) != len(np_type.names): + return value + return tuple( + _as_record_parts(np_type.fields[name][0], part) + for name, part in zip(np_type.names, value, strict=True) + ) def _leaf_arrays(arr): diff --git a/graphblas/core/utils.py b/graphblas/core/utils.py index 74d6c8344..48da77a5f 100644 --- a/graphblas/core/utils.py +++ b/graphblas/core/utils.py @@ -293,7 +293,24 @@ def normalize_values(self, values, dtype, shape=None, is_iso=False): if dtype is not None: dtype = lookup_dtype(dtype) if dtype != self.dtype: - values = values.astype(dtype.np_type) # copies + np_type = dtype.np_type + element_shape = self.dtype.np_type.shape + if element_shape and np_type.subdtype is not None: + # ``values`` has each element's own axes last already. astype to the + # array dtype would add them again, repeating each number into one. + from ..exceptions import DomainMismatch + from .udt import _strip_leading_ones + + if _strip_leading_ones(element_shape) != _strip_leading_ones(np_type.shape): + # As a store says it (``to_dense(fill_value, dtype)`` stores) + raise DomainMismatch( + f"cannot convert values of type {self.dtype} to {dtype}: array UDTs " + "must have the same shape, apart from leading axes of length 1" + ) + index_shape = values.shape[: values.ndim - len(element_shape)] + values = values.astype(np_type.base).reshape(index_shape + np_type.shape) + else: + values = values.astype(np_type) # copies else: dtype = self.dtype if is_iso: diff --git a/graphblas/core/vector.py b/graphblas/core/vector.py index 451be6ccb..17676ac51 100644 --- a/graphblas/core/vector.py +++ b/graphblas/core/vector.py @@ -2,7 +2,7 @@ import numpy as np -from .. import backend, binary, monoid, select, semiring, unary +from .. import backend, binary, indexunary, monoid, select, semiring, unary from ..dtypes import _INDEX, FP64, INT64, UINT64, lookup_dtype, unify from ..exceptions import DimensionMismatch, GrB_NO_VALUE, NoValue, check_status, check_status_carg from . import _supports_udfs, automethods, ffi, lib, utils @@ -111,6 +111,23 @@ def _select_mask(updater, obj, mask): updater << obj.dup(mask=mask) +def _udt_select_comparison(op, dtype): + """``binary.eq`` or ``ne`` for select op ``valueeq`` or ``valuene`` on UDT ``dtype``, else None. + + No select op lifts to UDTs, while ``==`` and ``!=`` do, so such a select keeps + the values they find equal or unequal to the thunk (:func:`_select_by_comparison`). + """ + if dtype._is_udt: + for name, comparison in [("valueeq", binary.eq), ("valuene", binary.ne)]: + if op is getattr(select, name) or op is getattr(indexunary, name): + return comparison + return None + + +def _select_by_comparison(updater, obj, mask_expr): + _select_mask(updater, obj, mask_expr.new(name="M_temp").V) + + def _abs_type(dtype): return dtype if dtype.np_type.kind == "c" else FP64 @@ -1738,9 +1755,10 @@ def select(self, op, thunk=None): if thunk is None: thunk = False # most basic form of 0 when unifying dtypes given = thunk + comparison = _udt_select_comparison(op, self.dtype) if type(thunk) is not Scalar: # A thunk beside a UDT becomes an element of it, so it must fit, as in apply. - thunk, dtype = _literal_operand(self.dtype, thunk, op) + thunk, dtype = _literal_operand(self.dtype, thunk, comparison or op) try: thunk = Scalar.from_value(thunk, dtype, is_cscalar=None, name="") except TypeError: @@ -1752,6 +1770,16 @@ def select(self, op, thunk=None): extra_message="Literal scalars also accepted.", op=op, ) + if comparison is not None: + mask_expr = self.apply(comparison, right=thunk) # typed now, so errors raise now + return VectorExpression( + method_name, + None, + [self, thunk, _select_by_comparison, (self, mask_expr)], + expr_repr=f"{{0.name}}.select({op}, thunk={{1._expr_name}})", + size=self._size, + dtype=self.dtype, + ) op = get_typed_op(op, self.dtype, thunk.dtype, is_right_scalar=True, kind="select") op, thunk = _value_select_in_own_type(op, self.dtype, thunk, given) self._expect_op(op, ("SelectOp", "IndexUnaryOp"), within=method_name, argname="op") diff --git a/graphblas/select/__init__.py b/graphblas/select/__init__.py index 43e8c206e..ed59d805b 100644 --- a/graphblas/select/__init__.py +++ b/graphblas/select/__init__.py @@ -45,18 +45,25 @@ def _match_expr(parent, expr): The argument must match the parent, so this _won't_ be rewritten: ``A.select(B < 5)`` """ + from ..core.utils import _Pointer + args = expr.args op = expr.op if ( expr.method_name == "apply" and len(args) == 2 and op.parent in _binary_to_select - and args[1]._is_scalar and args[0] is parent - and op.type in _binary_to_select[op.parent] ): - selectop = _binary_to_select[op.parent][op.type] - return selectop(parent, args[1]) + thunk = args[1] + if op.type._is_udt and op.parent.name in {"eq", "ne"}: + # ``select`` keeps what ``==`` or ``!=`` finds on a UDT, whose Scalar may + # be behind a pointer. + thunk = thunk.val if isinstance(thunk, _Pointer) else thunk + return _binary_to_select[op.parent](parent, thunk) + if thunk._is_scalar and op.type in _binary_to_select[op.parent]: + selectop = _binary_to_select[op.parent][op.type] + return selectop(parent, thunk) def value(expr): diff --git a/graphblas/tests/test_dtype.py b/graphblas/tests/test_dtype.py index 5b2a565b0..aca03cad5 100644 --- a/graphblas/tests/test_dtype.py +++ b/graphblas/tests/test_dtype.py @@ -375,8 +375,10 @@ def test_register_anonymous_zero_size(): for dtype in [str, bytes, np.dtype("U"), np.dtype([]), np.dtype((np.int64, (0,)))]: with pytest.raises(ValueError, match="size of 0"): dtypes.register_anonymous(dtype) - with pytest.raises(ValueError, match="Unknown dtype"): + with pytest.raises(ValueError, match="Unknown dtype") as exc_info: dtypes.lookup_dtype(str) + # The reason it cannot be registered is chained; it was dropped. + assert "size of 0" in str(exc_info.value.__cause__) def test_register_anonymous_array_numba_cannot_type(): @@ -400,6 +402,28 @@ def test_register_anonymous_array_numba_cannot_type(): dtypes.register_anonymous(np_type) assert isinstance(exc_info.value.__cause__, AttributeError) assert np_type not in core.dtypes._registry + # So does a lookup, which registers a dtype, with that reason as its cause. + with pytest.raises(ValueError, match="Unknown dtype") as exc_info: + gb.Vector(np_type, size=1) + assert f"cannot type an array of {element}" in str(exc_info.value.__cause__) + + +def test_register_anonymous_raw_bytes_numba_cannot_type(): + # A void dtype with no fields is raw bytes, which Numba cannot type, alone or in + # an array or record. Its own error was "NumbaNotImplementedError: |V4" in color. + for np_type, raw in [ + (np.dtype("V4"), "V4"), + (np.dtype(("V4", (2,))), "V4"), + (np.dtype([("_rbt_a", "V3"), ("_rbt_b", np.int8)]), "V3"), + (np.dtype([("_rbt_c", [("_rbt_d", "V2")])]), "V2"), + ]: + if not supports_udfs: + assert dtypes.register_anonymous(np_type).np_type == np_type + continue + match = rf"cannot type raw bytes \(\|{raw}\).*np.dtype\(\(np.uint8, {raw[1:]}\)\)" + with pytest.raises(ValueError, match=match): + dtypes.register_anonymous(np_type) + assert np_type not in core.dtypes._registry def test_register_anonymous_from_dataclass_deferred_annotations(): diff --git a/graphblas/tests/test_indexbinary.py b/graphblas/tests/test_indexbinary.py index 196787d17..3b089f210 100644 --- a/graphblas/tests/test_indexbinary.py +++ b/graphblas/tests/test_indexbinary.py @@ -85,6 +85,45 @@ def add_theta(x, ix, jx, y, iy, jy, theta): # pragma: no cover (numba) assert list(C.to_coo()[2]) == [18, 19] +def test_int_theta_beyond_its_type(): + """An int theta its type cannot hold names the int and the type, as ``apply`` does. + + It was cffi's bare ``OverflowError: int too big to convert``. + """ + + def add_theta(x, ix, jx, y, iy, jy, theta): # pragma: no cover (numba) + return x + y + theta + + op = indexbinary.register_anonymous(add_theta) + for theta in [2**63, -(2**63) - 1, 2**70]: + with pytest.raises( + OverflowError, match=f"^Python integer {theta} out of bounds for int64," + ): + op(theta) + with pytest.raises(OverflowError, match=f"^Python integer {2**64} out of bounds for uint64$"): + op(2**64, dtype="UINT64") + v = Vector.from_coo([0], [1], dtype="UINT64") + assert v.ewise_mult(v, op(2**63, dtype="UINT64")).new()[0].new() == 2**63 + 2 + assert v.ewise_mult(v, op(2**70, dtype="FP64")).new()[0].new() == 2.0**70 + + +def test_index_with_too_many_types(): + """``op[a, b, c, d]`` names what the op takes; it was ValueError from unpacking a pair.""" + + def add_theta(x, ix, jx, y, iy, jy, theta): # pragma: no cover (numba) + return x + y + theta + + op = indexbinary.register_anonymous(add_theta) + types = ("INT64",) * 4 + for parent in [op, op(1)]: + with pytest.raises( + TypeError, match=r"takes a dtype, or a tuple of two \(x and y\) .* of 4" + ): + parent[types] + assert types not in op + assert types[:3] in op + + def test_untyped_call_is_typed_when_it_runs(): """``op(theta)`` takes its type from both inputs and theta each time it runs. diff --git a/graphblas/tests/test_matrix.py b/graphblas/tests/test_matrix.py index 91235a307..bd981b8d7 100644 --- a/graphblas/tests/test_matrix.py +++ b/graphblas/tests/test_matrix.py @@ -4225,7 +4225,7 @@ def assert_same(A, B): A[1, 2] = (1, 5.6) assert_same(A, Matrix.ss.deserialize(A.ss.serialize(), dtype=long_udt)) # An anonymous one needs dtype= (see test_vector.py's test_udt for the errors) - with pytest.raises(SyntaxError if hasattr(lib, "GxB_Serialized_get_String") else ValueError): + with pytest.raises(ValueError, match="pass it as dtype="): Matrix.ss.deserialize(A.ss.serialize()) # A named one is found by its name named_long_dtype = np.dtype([("x", np.bool_), ("z" * 1000, np.float64)], align=False) @@ -4837,6 +4837,7 @@ def test_subarray_dtypes(): assert A.isequal(B, check_dtype=True) for method in ["coo", "csr", "csc", "dcsr", "dcsc", "edgelist"]: B = getattr(A, f"from_{method}")(*getattr(A, f"to_{method}")()) + assert A.isequal(B, check_dtype=True) B = Matrix.from_dicts(A.to_dicts()) assert A.isequal(B, check_dtype=True) B = Matrix.from_dicts(A.to_dicts(), A.dtype) @@ -4859,12 +4860,170 @@ def test_subarray_dtypes(): Full2 = Matrix.from_dense(b2) assert Full1.isequal(Full2, check_dtype=True) Full2 = Matrix.from_dense(Full1.to_dense()) - assert Full2.isequal(Full2, check_dtype=True) + assert Full1.isequal(Full2, check_dtype=True) if suitesparse: Full2 = Matrix.ss.import_fullr(b2) assert Full1.isequal(Full2, check_dtype=True) +def _udt_entries(M): + rows, cols, values = M.to_coo() + order = np.lexsort((cols, rows)) + return rows[order], cols[order], values[order] + + +def _assert_same_udt_entries(A, B): + # isequal needs Numba for a UDT, so compare the entries with numpy + assert B.dtype == A.dtype + assert B.shape == A.shape + for a, b in zip(_udt_entries(A), _udt_entries(B), strict=True): + np.testing.assert_array_equal(b, a) + + +@pytest.mark.skipif("not suitesparse") +@pytest.mark.parametrize("shape", [(4,), (2, 2)]) +def test_subarray_dtype_ss_column_order(shape): + # values[i, j] is the element at (i, j) in every dense format, exported and + # imported. Column order used to put the element's own axes slowest, both ways, + # so the numbers came out scrambled and went in scrambled, and an unpack then + # pack round trip only worked because the two cancelled. + expected = np.arange(2 * 3 * 4, dtype=np.int64).reshape(2, 3, *shape) + A = Matrix.from_dense(expected) + assert A.dtype.np_type.shape == shape + bitmap = np.ones((2, 3), dtype=bool) + for fmt in ["fullr", "fullc", "bitmapr", "bitmapc"]: + info = A.ss.export(fmt) + np.testing.assert_array_equal(info["values"], expected) + _assert_same_udt_entries(A, Matrix.ss.import_any(**info)) + kwargs = {"bitmap": bitmap} if fmt.startswith("bitmap") else {} + import_ = getattr(Matrix.ss, f"import_{fmt}") + for values in [expected, np.asfortranarray(expected), expected.tolist()]: + _assert_same_udt_entries(A, import_(values=values, **kwargs)) + B = import_(values=values, dtype=A.dtype, **kwargs) + _assert_same_udt_entries(A, B) + values = expected.copy() + kwargs = {key: val.copy() for key, val in kwargs.items()} # may be taken + B = import_(values=values, take_ownership=True, **kwargs) + _assert_same_udt_entries(A, B) + if fmt.endswith("c"): + # Never taken: the entries must be swapped into a copy + assert values.flags.owndata + assert values.flags.writeable + B = A.dup() + getattr(B.ss, f"pack_{fmt}")(**B.ss.unpack(fmt)) + _assert_same_udt_entries(A, B) + # Numbers given flat, in column order, still import as they are + flat = expected.swapaxes(0, 1).ravel() + B = Matrix.ss.import_fullc(values=flat, nrows=2, ncols=3, dtype=A.dtype) + _assert_same_udt_entries(A, B) + raw = A.ss.export("fullc", raw=True) # may have spare values after the six + np.testing.assert_array_equal(raw["values"][:6].ravel(), flat) + # A list of elements in column order, as ``raw=True`` exports them, imports as it is + for fmt in ["fullc", "bitmapc"]: + info = A.ss.export(fmt, raw=True) + assert info["values"].shape[1:] == shape + _assert_same_udt_entries(A, Matrix.ss.import_any(**info)) + B = A.dup() + B.ss.pack_any(**B.ss.unpack(fmt, raw=True)) + _assert_same_udt_entries(A, B) + kwargs = {"bitmap": np.ones(6, dtype=bool)} if fmt == "bitmapc" else {} + values = expected.swapaxes(0, 1).reshape(6, *shape) + B = getattr(Matrix.ss, f"import_{fmt}")( + values=values, nrows=2, ncols=3, dtype=A.dtype, take_ownership=True, **kwargs + ) + _assert_same_udt_entries(A, B) + assert not values.flags.owndata # taken, as it is + # Element axes that fit neither layout are refused, not read as a list of elements + # (for INT64[2, 2], values of shape (2, 3, 4) put [8, 9, 10, 11] at (0, 1)) + bad = expected.reshape(2, 3, 4) if len(shape) > 1 else expected.reshape(2, 12) + kwargs = {"bitmap": np.ones((2, 3), dtype=bool)} if fmt == "bitmapc" else {} + with pytest.raises(ValueError, match=r"in column order must have shape \(nrows, ncols"): + getattr(Matrix.ss, f"import_{fmt}")( + values=bad, nrows=2, ncols=3, dtype=A.dtype, **kwargs + ) + # An iso value given as one element (for INT64[2, 2], it was read transposed) + B = getattr(Matrix.ss, f"import_{fmt}")( + values=expected[0, 1], nrows=2, ncols=3, dtype=A.dtype, is_iso=True, **kwargs + ) + np.testing.assert_array_equal(B.to_dense(), np.broadcast_to(expected[0, 1], (2, 3, *shape))) + # Elements missing from a bitmap are left out + del A[0, 1] + info = A.ss.export("bitmapc") + mask = info["bitmap"] + assert mask.sum() == A.nvals == 5 + np.testing.assert_array_equal(info["values"][mask], expected[mask]) + _assert_same_udt_entries(A, Matrix.ss.import_any(**info)) + + +@pytest.mark.skipif("not suitesparse") +def test_subarray_dtype_ss_export_spare_values(): + # A dense export of an array UDT given more values than entries raised + # ValueError, and unpack left the Matrix empty. SuiteSparse 7 always gives more + # (its memory pool rounds allocations up); here a full Matrix resized smaller + # keeps its larger buffer, in either layout. + big = np.arange(3 * 3 * 4, dtype=np.int64).reshape(3, 3, 4) + for layout, nrows, ncols in [("by_row", 2, 3), ("by_col", 3, 2)]: + expected = big[:nrows, :ncols] + + def resized(layout=layout, nrows=nrows, ncols=ncols): + A = Matrix.from_dense(big) + A.ss.config["format"] = layout + A.resize(nrows, ncols) + return A + + for fmt in ["fullr", "fullc", "bitmapr", "bitmapc"]: + A = resized() + np.testing.assert_array_equal(A.to_dense(), expected) + info = A.ss.unpack(fmt) + np.testing.assert_array_equal(info["values"], expected) + getattr(A.ss, f"pack_{fmt}")(**info) + np.testing.assert_array_equal(A.to_dense(), expected) + info = resized().ss.export(fmt, give_ownership=True) + np.testing.assert_array_equal(info["values"], expected) + np.testing.assert_array_equal(Matrix.ss.import_any(**info).to_dense(), expected) + + +@pytest.mark.skipif("not suitesparse") +def test_subarray_dtype_ss_row_order_round_trip_takes_the_buffer(): + # Unpacked in row order from a buffer of exactly its size (256 bytes, which + # SuiteSparse 7's pool does not round up), an array UDT's values own that buffer, + # so packing them back with take_ownership=True hands it over without a copy + expected = np.arange(2 * 4 * 4, dtype=np.int64).reshape(2, 4, 4) + for fmt in ["fullr", "bitmapr"]: + A = Matrix.from_dense(expected) + info = A.ss.unpack(fmt) + values = info["values"] + assert values.flags.owndata + getattr(A.ss, f"pack_{fmt}")(take_ownership=True, **info) + assert not values.flags.owndata + np.testing.assert_array_equal(A.to_dense(), expected) + + +# Without SuiteSparse, building a UDT uses dup_op=binary.any, which needs Numba. +@pytest.mark.skipif("not suitesparse and not supports_udfs") +def test_subarray_dtype_values_in_another_array_dtype(): + # Asked for another array dtype of the same shape, the values convert number by + # number. astype to the array dtype repeated each number into an element: shape + # (n, 4, 4), for an iso Matrix too. + dense = np.arange(2 * 3 * 4, dtype=np.int64).reshape(2, 3, 4) + A = Matrix.from_dense(dense) + iso = Matrix("INT64[4]", nrows=2, ncols=3) + iso[:, :] = [0, 1, 2, 3] + for M, expected in [(A, dense), (iso, np.broadcast_to([0, 1, 2, 3], (2, 3, 4)))]: + got = M.to_dense(dtype="FP64[4]") + assert got.dtype == np.float64 + np.testing.assert_array_equal(got, expected) + np.testing.assert_array_equal(M.to_dense(dtype="FP32[1, 4]"), expected[:, :, None, :]) + with pytest.raises(DomainMismatch, match="must have the same shape"): + M.to_dense(dtype="FP64[2]") + values = expected.reshape(6, 4) + for method in ["coo", "csr", "dcsr"]: + np.testing.assert_array_equal(getattr(M, f"to_{method}")(dtype="FP64[4]")[-1], values) + values = expected.swapaxes(0, 1).reshape(6, 4) + for method in ["csc", "dcsc"]: + np.testing.assert_array_equal(getattr(M, f"to_{method}")(dtype="FP64[4]")[-1], values) + + # Without SuiteSparse, building a UDT uses dup_op=binary.any, which needs Numba. @pytest.mark.skipif("not suitesparse and not supports_udfs") def test_subarray_dtype_values_shape(): @@ -4938,6 +5097,29 @@ def entries(M): assert entries(B) == entries(Matrix.from_coo([0, 2, 1], cols.ravel(), [1, 2, 3])) +@pytest.mark.skipif("not suitesparse") +def test_ss_export_coo_of_an_iso_udt(): + # An iso Matrix exports as "coo" with one value. For a UDT this called + # GrB_Matrix_extractTuples_, which does not exist (AttributeError), and the + # value came from a reduce, which a UDT can do only with Numba. + arr = dtypes.register_anonymous(np.dtype((np.int64, (4,))), "_IsoCooArr") + rec = dtypes.register_anonymous(np.dtype([("a", np.int64), ("b", np.float64)]), "_IsoCooRec") + A = Matrix.ss.import_fullr(values=[[1, 2, 3, 4]], nrows=2, ncols=3, dtype=arr, is_iso=True) + R = Matrix.ss.import_fullr( + values=np.array([(1, 2.5)], rec.np_type), nrows=2, ncols=3, is_iso=True + ) + F = Matrix.ss.import_fullr(values=[2.5], nrows=2, ncols=3, is_iso=True) + for M, value in [(A, [[1, 2, 3, 4]]), (R, [(1, 2.5)]), (F, [2.5])]: + assert M.ss.is_iso + info = M.ss.export("coo") + assert info["is_iso"] + assert info["values"].tolist() == value # one value, (1, 4) for the array UDT + B = Matrix.ss.import_any(**info) + assert B.dtype == M.dtype + for got, want in zip(B.to_coo(), M.to_coo(), strict=True): + np.testing.assert_array_equal(got, want) + + @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 @@ -4972,7 +5154,8 @@ def test_ss_import_values_size(): ) with pytest.raises(ValueError, match=match): Matrix.from_dcsr([0, 1], [0, 1, 2], [0, 1], np.ones(2), "FP64[2]") - with pytest.raises(ValueError, match=r"have 4 or more elements of dtype \S+, not 2"): + match = r"values of \S+ in column order must have shape \(nrows, ncols, 2\), an element" + with pytest.raises(ValueError, match=match): Matrix.ss.import_fullc(values=np.ones((2, 2, 1)), dtype="FP64[2]") # A scalar is not a whole element of a sub-array dtype (it said "not 0") with pytest.raises(ValueError, match=r"a scalar does not fill one element of dtype \S+, of"): diff --git a/graphblas/tests/test_op.py b/graphblas/tests/test_op.py index 57c14c215..94d1ae78d 100644 --- a/graphblas/tests/test_op.py +++ b/graphblas/tests/test_op.py @@ -169,6 +169,12 @@ def test_get_typed_op(): with pytest.raises(AttributeError): # get_typed_op expects dtypes to already be dtypes operator.get_typed_op(binary.plus, dtypes.INT64, "bad dtype") + # A tuple of types other than a pair (it was ValueError from unpacking it) + for parent in [binary.plus, monoid.plus, indexunary.rowle]: + for types in [("INT64",) * 3, ("INT64",)]: + with pytest.raises(TypeError, match=rf"takes a dtype or a pair .* of {len(types)}$"): + parent[types] + assert types not in parent @pytest.mark.skipif("supports_udfs") @@ -213,6 +219,30 @@ def test_unaryop_on_udt_needs_numba(): v.apply(unop).new() +@pytest.mark.skipif("supports_udfs") +def test_udt_isequal_needs_numba(): + """Without Numba, ``isequal`` of a Vector or Matrix of a UDT says it needs Numba. + + It raised binary.eq's bare KeyError. A Scalar compares in numpy. + """ + rec = dtypes.register_anonymous(np.dtype([("nn_c", np.int8), ("nn_d", np.float64)])) + text = dtypes.register_anonymous(np.dtype("S4"), "_NoNumbaText") + v = Vector(rec, size=2) + v[0] = (1, 2.0) + t = Matrix(text, nrows=1, ncols=1) + t[0, 0] = b"ab" + for expr in [ + lambda: v.isequal(v), + lambda: t.isequal(t), + lambda: v == v, + lambda: v.select("==", (1, 2.0)), + ]: + with pytest.raises(KeyError, match="work on UDTs only with Numba installed"): + expr() + assert gb.Scalar.from_value((1, 2.0), rec).isequal(v[0].new()) + assert gb.Scalar.from_value(b"ab", text) == b"ab" + + @pytest.mark.skipif("not supports_udfs") def test_unaryop_udf(): def plus_one(x): @@ -2161,6 +2191,26 @@ def _sum_a(x, y): assert record_udt not in sum_a_monoid +@pytest.mark.skipif("not supports_udfs") +def test_udt_monoid_on_a_builtin_type(): + """A monoid over an ``is_udt=True`` op works on a built-in type the op returns. + + GrB_Monoid_new_UDT refused the op typed FP64 x FP64 -> FP64 with DomainMismatch + and no message; the typed GrB_Monoid_new_FP64 takes it. + """ + op = BinaryOp.register_anonymous(lambda x, y: x + y, is_udt=True) + mon = Monoid.register_anonymous(op, np.float64(0.0)) + for dtype in [FP64, FP32]: + assert mon[dtype].return_type == dtype + assert Vector.from_coo([0, 1, 2], [1, 2, 3], dtype=dtype).reduce(mon).new() == 6 + assert Vector(dtype, 3).reduce(mon, allow_empty=False).new() == 0 + A = Matrix.from_coo([0, 0, 1], [0, 1, 1], [1.0, 2.0, 3.0]) + assert A.reduce_rowwise(mon).new().isequal(Vector.from_coo([0, 1], [3.0, 3.0])) + # An op that returns another type is still refused up front + with pytest.raises(TypeError, match="must return its input type, INT8, not INT64"): + mon[INT8] + + @pytest.mark.skipif("not supports_udfs") @pytest.mark.slow def test_udt_tuple_return_semiring(record_udt): @@ -2262,6 +2312,56 @@ def _array(x, y): # pragma: no cover (numba; raises before execution) BinaryOp.register_anonymous(func, is_udt=True)[udt] +@pytest.mark.skipif("not supports_udfs") +def test_udt_udf_reading_a_missing_field(): + """``x.name`` for a field the record lacks names the op and types, as ``x["name"]`` does. + + Numba types the attribute form with a bare ``KeyError: 'name'``; before 0.60, + with a TypingError that holds it, and 0.58 and 0.59 also warn about that. + ``x["name"]``'s "Field" line comes quoted on 0.60 (in a bare KeyError) and before. + """ + rec = dtypes.register_anonymous(np.dtype([("mf_a", "i1"), ("mf_b", "f8")]), "_MissFieldRec") + other = dtypes.register_anonymous(np.dtype([("mf_c", "i1"), ("mf_d", "f8")]), "_MissFieldRec2") + unknown = BinaryOp.register_anonymous(lambda x, y: x.mf_zzz, "_miss_field", is_udt=True) + with pytest.raises( + UdfParseError, + match=r"_miss_field does not work with \(_MissFieldRec, _MissFieldRec\): " + r"Field 'mf_zzz' was not found", + ): + unknown[rec] + assert rec not in unknown + mixed = BinaryOp.register_anonymous(lambda x, y: (x.mf_a, y.mf_d), "_miss_mix", is_udt=True) + with pytest.raises(UdfParseError, match=r"\(_MissFieldRec2, _MissFieldRec\): Field 'mf_a'"): + mixed[other, rec] + unary_op = UnaryOp.register_anonymous(lambda x: x.mf_zzz, "_miss_field_u", is_udt=True) + with pytest.raises(UdfParseError, match=r"unary._miss_field_u does not work with"): + unary_op[rec] + by_key = UnaryOp.register_anonymous(lambda x: x["mf_zzz"], "_miss_field_k", is_udt=True) + with pytest.raises( + UdfParseError, + match=r"unary._miss_field_k does not work with \(_MissFieldRec\): " + r"Field 'mf_zzz' was not found in record with fields \('mf_a', 'mf_b'\)$", + ): + by_key[rec] + # The TypingError from Numba before 0.60, as 0.59.1 words it, whatever Numba runs here + from numba import TypingError + + from graphblas.core.operator.udf import _missing_field, _summarize_numba_typing_error + + old_style = ( + "Failed in nopython mode pipeline (step: nopython frontend)\n" + 'Internal error at resolving type of attribute "mf_zzz" of "x".\n' + "\x1b[1m\x1b[1m'{}'\x1b[0m\n\x1b[0m\x1b[1mDuring: typing of get attribute at f.py (1)" + ) + assert _missing_field(TypingError(old_style.format("mf_zzz"))) == "mf_zzz" + assert _missing_field(TypingError(old_style.format("other error"))) is None + assert _missing_field(KeyError("mf_zzz")) == "mf_zzz" + # Numba 0.60 raises x["mf_zzz"]'s sentence as a bare KeyError, which is no name + by_key_060 = KeyError("Field 'mf_zzz' was not found in record with fields ('mf_a', 'mf_b')") + assert _missing_field(by_key_060) is None + assert _summarize_numba_typing_error(by_key_060) == by_key_060.args[0] + + @pytest.mark.skipif("not supports_udfs") @pytest.mark.slow def test_udt_return_type_errors_array_udt(): @@ -2464,6 +2564,41 @@ def test_udt_unary_ainv_abs_array(array_udt): np.testing.assert_array_equal(abs_c[1].new().value, [4.0, 5.0, 6.0]) +@pytest.mark.skipif("not supports_udfs") +def test_udt_abs_of_complex_leaves_is_real(): + """``unary.abs`` types a complex leaf as on FC64 or FC32: FP64 or FP32, as numpy's abs. + + It kept the complex type, as ``(5+0j)``. + """ + c3 = dtypes.register_anonymous(np.dtype((np.complex128, (3,))), "_AbsC16x3") + c2 = dtypes.register_anonymous(np.dtype((np.complex64, (2,))), "_AbsC8x2") + rec = dtypes.register_anonymous( + np.dtype([("abs_c", np.complex128), ("abs_f", np.complex64), ("abs_i", np.int8)]), + "_AbsCRec", + ) + v = Vector(c3, size=1) + v[0] = [3 + 4j, -1, -2j] + w = Vector(c2, size=1) + w[0] = [3 + 4j, 1j] + r = Vector(rec, size=1) + r[0] = (3 + 4j, 6 - 8j, -3) + for x, expected_type in [ + (v, np.dtype((np.float64, (3,)))), + (w, np.dtype((np.float32, (2,)))), + (r, np.dtype([("abs_c", np.float64), ("abs_f", np.float32), ("abs_i", np.int8)])), + ]: + assert unary.abs[x.dtype].return_type.np_type == expected_type + assert x.apply(unary.abs).new().dtype.np_type == expected_type + np.testing.assert_array_equal(v.apply(unary.abs).new()[0].new().value, [5, 1, 2]) + np.testing.assert_array_equal(w.apply(unary.abs).new()[0].new().value, [5, 1]) + assert r.apply(unary.abs).new()[0].new().value.tolist() == (5.0, 10.0, 3) + assert unary.ainv[c3].return_type == c3 + # Stored into the complex type, as an FP64 result is + stored = Vector(c3, size=1) + stored << v.apply(unary.abs) + np.testing.assert_array_equal(stored[0].new().value, [5, 1, 2]) + + @pytest.mark.skipif("not supports_udfs") @pytest.mark.slow def test_udt_jit_typedef(): @@ -2686,6 +2821,35 @@ def _fieldwise_max(x, y): assert v.reduce(Monoid.register_anonymous(fieldwise_max, (lowest, lowest))).new() == (5, 2) +@pytest.mark.skipif("not supports_udfs") +def test_monoid_refuses_an_empty_scalar_identity(): + """An empty Scalar identity raises EmptyObject, as an empty theta does. + + A UDT monoid took it as zero bytes, so a field-wise product reduced to all + zeros, and a built-in type raised ``TypeError`` from converting None. + """ + udt = dtypes.register_anonymous(np.dtype([("ei_a", np.int8), ("ei_b", np.float64)])) + + def _fieldwise_times(x, y): # pragma: no cover (numba) + return (x["ei_a"] * y["ei_a"], x["ei_b"] * y["ei_b"]) + + times_udt = BinaryOp.register_anonymous(_fieldwise_times, "_ei_times", is_udt=True) + times_int = BinaryOp.register_anonymous(lambda x, y: x * y, "_ei_times_int") + for binop, identity in [ + (times_udt, gb.Scalar(udt)), + (times_udt, gb.Scalar(udt, is_cscalar=True)), + (times_int, gb.Scalar(INT64)), + (times_int, {INT64: gb.Scalar(INT64)}), + ]: + with pytest.raises(gb.exceptions.EmptyObject, match="identity is an empty Scalar"): + Monoid.register_anonymous(binop, identity) + v = Vector(udt, size=2) + v[0] = (2, 3.0) + v[1] = (5, 7.0) + one = gb.Scalar.from_value((1, 1.0), udt) + assert v.reduce(Monoid.register_anonymous(times_udt, one)).new() == (10, 21.0) + + @pytest.mark.skipif("not supports_udfs") @pytest.mark.parametrize("shape", [(4,), (3, 2)]) def test_udt_array_wrapper_stays_within_element(shape): @@ -5601,6 +5765,139 @@ def test_udt_array_literal_broadcasts_with_an_array_udt(): # A literal that must become an element must still fit it. with pytest.raises(ValueError, match="does not fit _BcastI8x3"): v.ewise_union(v, binary.plus, rows, 0) + # Shapes that do not broadcast raise the KeyError of two such UDTs, which + # names both; they raised numpy's "could not broadcast" from converting. + pair = np.array([1, 2]) + both = r"broadcast together; got (int8\[3\] vs int\d+\[2\]|int\d+\[2\] vs int8\[3\])" + for expr in [ + lambda: v + pair, + lambda: v == pair, + lambda: binary.plus(v, [1, 2]), + lambda: v.apply(binary.plus, left=pair), + lambda: v[0].new() == pair, + lambda: v.apply(monoid.plus, right=pair), + ]: + with pytest.raises(KeyError, match=both): + expr() + with pytest.raises(ValueError, match="does not fit _BcastI8x3"): + v.ewise_union(v, binary.plus, pair, 0) + # An empty array has no UDT layout; numpy says why it does not fit + for empty in [np.array([]), np.zeros((2, 0)), np.zeros((0, 3))]: + with pytest.raises(ValueError, match="could not broadcast"): + binary.plus(v, empty) + + +@pytest.mark.skipif("not supports_udfs") +def test_udt_record_literal_may_be_a_list(): + """A list laid out like a record converts as a tuple does, as udt.rst promises. + + numpy takes only a tuple for a record, so a list raised "could not broadcast + input array from shape (2,) into shape (1,)" from every op and store. + """ + rec = dtypes.register_anonymous(np.dtype([("rl_x", np.int8), ("rl_y", np.int8)]), "_ListRec") + nested = dtypes.register_anonymous( + np.dtype([("rl_p", [("rl_a", np.int8), ("rl_b", np.float64)]), ("rl_v", np.float64, (2,))]), + "_ListNested", + ) + v = Vector(rec, size=2) + v[0] = [1, 2] + v[1] = (3, 4) + assert v[0].new().value.tolist() == (1, 2) + assert gb.Scalar.from_value([1, 2], rec) == (1, 2) + assert binary.plus(v, [1, 2]).new()[0].new().value.tolist() == (2, 4) + float_y = np.dtype([("rl_x", np.int8), ("rl_y", np.float64)]) + assert binary.plus(v, [1, 2.5]).new().dtype.np_type == float_y + assert (v == [1, 2]).new().to_coo()[1].tolist() == [True, False] + assert gb.Scalar.from_value((1, 2), rec).isequal([1, 2]) + first = BinaryOp.register_anonymous(lambda x, y: y, "_list_rec_second", is_udt=True) + assert v.apply(first, right=[5, 6]).new()[1].new().value.tolist() == (5, 6) + # In a nested record too, where a list for an array field stays the field's value + n = Vector(nested, size=1) + n[0] = [[1, 2.0], [3.0, 4.0]] + for value, expected in [ + (n[0].new().value, [(1, 2.0), [3.0, 4.0]]), + (binary.plus(n, [[1, 1], [1, 1]]).new()[0].new().value, [(2, 3.0), [4.0, 5.0]]), + ]: + assert [value["rl_p"].tolist(), value["rl_v"].tolist()] == expected + assert (n == {"rl_p": [1, 2.0], "rl_v": [3.0, 4.0]}).new()[0].new().value + # A list of another length is refused by its length, as a tuple is + for expr in [lambda: binary.plus(v, [1, 2, 3]), lambda: gb.Scalar.from_value([1, 2, 3], rec)]: + with pytest.raises(ValueError, match="tuple of length 3 to structure with 2 fields"): + expr() + + +@pytest.mark.skipif("not supports_udfs") +def test_udt_select_by_equality(): + """``select("==", x)`` and ``"!="`` keep what ``==`` and ``!=`` find, as for built-in dtypes. + + No select op lifts to UDTs, so they raised ``valueeq does not work with``, + and ``v.select(v == x)`` raised AttributeError from a UDT Scalar's pointer. + """ + rec = dtypes.register_anonymous(np.dtype([("se_a", np.int8), ("se_b", np.float64)]), "_SelEq") + v = Vector(rec, size=4) + v[0] = (1, 2.0) + v[1] = (3, 4.0) + v[2] = (1, 2.0) + x = (1, 2.0) + for expr, indices in [ + (lambda: v.select("==", x), [0, 2]), + (lambda: v.select("!=", x), [1]), + (lambda: v.select(select.valueeq, gb.Scalar.from_value((3, 4.0), rec)), [1]), + (lambda: v.select(indexunary.valueeq, x), [0, 2]), + (lambda: select.valueeq(v, x), [0, 2]), + (lambda: select.value(v == x), [0, 2]), + (lambda: v.select(v == x), [0, 2]), + # Typed as == types it: an int beyond a field, or a float, equals nothing + (lambda: v.select("==", (300, 2.0)), []), + (lambda: v.select("==", (1.5, 2.0)), []), + ]: + result = expr().new() + assert result.dtype == rec + assert result.to_coo()[0].tolist() == indices + assert v.select("==", x).new()[0].new().value.tolist() == x + w = Vector(rec, size=4) + w[0] = (1, 1.0) + w(accum=binary.plus) << v.select("==", x) + assert w.to_coo()[1].tolist() == [(2, 3.0), (1, 2.0)] + # On a matrix that is not square, with entries off the diagonal, so that a lost + # transpose would show + A = Matrix(rec, nrows=2, ncols=3) + A[0, 2] = x + A[1, 0] = (5, 6.0) + for expr, expected in [ + (lambda: A.select("==", x), ([0], [2])), + (lambda: A.select("!=", x), ([1], [0])), + (lambda: A.select(A == x), ([0], [2])), + (lambda: A.T.select("==", (5, 6.0)), ([0], [1])), + ]: + rows, cols, _ = expr().new().to_coo() + assert (rows.tolist(), cols.tolist()) == expected + # A.select(A.T == x) was rewritten to A.select("==", x); it selects by the mask now + with gb.config.set(autocompute=True): + rows, cols, _ = A.T.select(A.T == x).new().to_coo() + assert (rows.tolist(), cols.tolist()) == ([2], [0]) + with pytest.raises(DimensionMismatch): + A.select(A.T == x).new() + # A mask and replace act on the selected entries as for any select + m = Matrix(rec, nrows=2, ncols=3) + m[0, 2] = (9, 9.0) + m[1, 1] = (9, 9.0) + m(Matrix.from_coo([0, 1], [2, 1], [True, True], nrows=2, ncols=3).S, replace=True) << ( + A.select("==", x) + ) + rows, cols, _ = m.to_coo() + assert (rows.tolist(), cols.tolist()) == ([0], [2]) + assert m[0, 2].new().value.tolist() == x + # Ordering does not lift to UDTs + with pytest.raises(KeyError, match="valuelt does not work with _SelEq"): + v.select("<", x) + # Errors raise at the call, as for any select, not when the expression runs + other = dtypes.register_anonymous(np.dtype([("se_c", np.int8), ("se_d", np.float64)])) + with pytest.raises(KeyError, match="record UDTs must share field names"): + v.select("==", gb.Scalar.from_value((1, 2.0), other)) + for thunk in [Vector.from_coo([0], [1]), object()]: + with pytest.raises(TypeError, match="Bad type for keyword argument `thunk=`"): + v.select("==", thunk) def test_udt_record_date_literals_compare_as_numpy(): @@ -5905,11 +6202,24 @@ def test_udt_store_casts_each_element_as_builtin_types_do(): C << A.mxm(A, semiring.plus_times) with pytest.raises(DomainMismatch, match="would need a temporary copy"): w << v.ewise_add(v, binary.plus) - # Only a built-in op that lifts to UDTs can, so not first or a user's op. + # Only a built-in op can, so not a user's op, nor any, which picks either operand. udf = BinaryOp.register_anonymous(lambda x, y: x, "_store_cast_first", is_udt=True) - for not_lifted in [binary.first, udf]: + for not_lifted in [binary.any, udf]: with pytest.raises(DomainMismatch, match="would need a temporary copy"): w << v.ewise_mult(v, not_lifted) + # identity, first and second store an operand, cast as a store casts it. + for u, expected in [(v, v.dup(i8s)), (x, z)]: + for expr in [ + u.apply(unary.identity), + u.ewise_mult(u, binary.first), + u.ewise_mult(u, binary.second), + u.apply(binary.first, right=u[0].new()), + u.apply(binary.second, left=u[0].new()), + ]: + stored = Vector(expected.dtype, size=u.size) + stored << expr + # As bytes, since z holds NaN + assert stored.to_coo()[1].tobytes() == expected.to_coo()[1].tobytes() with pytest.raises(DomainMismatch, match="extraction cannot"): w << v[[0, 1, 2]] # The computed result, stored as an existing object, casts in one pass. @@ -6382,6 +6692,19 @@ def test_udt_errors_name_the_real_problem(): expr().new() with pytest.raises(KeyError, match="lt does not work with _ErrI8"): v.apply(binary.lt, right=2.5) + # Also before a literal that does not fit an int8 field fails to convert. + rec = dtypes.register_anonymous(np.dtype([("err_a", "i1"), ("err_b", "f8")]), "_ErrRec") + r = Vector(rec, size=2) + s = gb.Scalar.from_value((1, 2.0), rec) + for expr in [ + lambda: v < 300, + lambda: r < (300, 2.0), + lambda: s < (300, 2.0), + lambda: r.apply(binary.lt, left=(300, 2.0)), + lambda: r.select("<", (300, 2.0)), + ]: + with pytest.raises(KeyError, match="lt does not work with _Err"): + expr() @pytest.mark.skipif("not supports_udfs") @@ -6424,11 +6747,11 @@ def test_udt_literal_rules_edge_cases(): with pytest.raises(ValueError, match=r"Fraction\(1, 2\) does not fit _EdgeI8"): v.apply(udf, right=Fraction(1, 2)) # A literal not laid out like an element (ragged, too long or too short for - # a field, not broadcastable), or not a number, has no type of its own, so - # it is converted into the UDT, and numpy says why it cannot be. A record - # says so itself for a part that does not broadcast to its array field. - # So too for a comparison, which first looks for an int beyond its leaf, a - # part holding numpy values, and text a bytes field cannot hold. + # a field), or not a number, has no type of its own, so it is converted + # into the UDT, and numpy says why it cannot be. A record says so itself + # for a part that does not broadcast to its array field. + # So too for a comparison, which first looks for an int beyond its leaf and + # a part holding numpy values. rec = dtypes.register_anonymous( np.dtype([("edge_v", np.float64, (3,)), ("edge_n", np.int8)]), "_EdgeRec" ) @@ -6441,7 +6764,6 @@ def test_udt_literal_rules_edge_cases(): text = dtypes.register_anonymous(np.dtype([("edge_s", "S3"), ("edge_k", np.int8)]), "_EdgeText") for expr, reason in [ (lambda: binary.plus(v, [1, [2, 3]]), "with a sequence"), - (lambda: v + np.array([1, 2]), "could not broadcast"), (lambda: v + "abc", "invalid literal"), (lambda: v.apply(udf, right="abc"), "invalid literal"), (lambda: binary.plus(r, (1, 2, 3)), "tuple of length 3"), @@ -6451,10 +6773,13 @@ def test_udt_literal_rules_edge_cases(): (lambda: binary.plus(r, ([np.float64(1), [2, 3]], 4)), "with a sequence"), (lambda: binary.plus(r, ([np.zeros((3, 2)), np.zeros((3, 1))], 4)), "with a sequence"), (lambda: binary.plus(r, (np.array(["a", "b", "c"]), 4)), "could not convert string"), - (lambda: Vector(text, size=1) == ("\xe9", 1), "can't encode"), ]: with pytest.raises(ValueError, match=reason): expr() + # Text a bytes field cannot hold equals nothing, as in numpy, rather than raise. + t = Vector(text, size=1) + t[0] = (b"\xe9", 1) + assert not (t == ("\xe9", 1)).new()[0].new().value c = Vector(c8s, size=1) c[0] = [1, 2, 3] assert c.apply(binary.plus, right=(1 + 2j, 0, 0)).new()[0].new().value[0] == 2 + 2j @@ -6873,6 +7198,20 @@ def test_udt_eq_compares_a_udt_that_is_neither_record_nor_array(): assert (four != b"abcdef").new()[0].new().value assert gb.Scalar.from_value(b"abcd", text) != b"abcdef" assert gb.Scalar.from_value(b"ab", text) == gb.Scalar.from_value(b"ab", wide) + # A str that bytes cannot hold (numpy encodes ASCII) equals no element, as in + # numpy, even one of its UTF-8 or Latin-1 bytes; it was UnicodeEncodeError. + latin = Vector(text, size=3) + latin[0] = b"\xe9" + latin[1] = "\xe9".encode() + for literal in ["\xe9", "a\xe9", np.str_("\xe9")]: + assert (latin == literal).new().to_coo()[1].tolist() == [False, False] + assert (latin != literal).new().to_coo()[1].tolist() == [True, True] + assert gb.Scalar.from_value(b"\xe9", text) != literal + assert not gb.Scalar.from_value(b"\xe9", text).isequal(literal) + text_rec = dtypes.register_anonymous(np.dtype([("eqo_s", "S4"), ("eqo_n", np.int8)])) + assert gb.Scalar.from_value((b"\xe9", 1), text_rec) != ("\xe9", 1) + with pytest.raises(ValueError, match="does not fit _EqOpaqueBytes"): + latin.apply(binary.first, right="\xe9") # A Scalar's == compares in numpy, so an element Numba has no type for compares too. half = dtypes.register_anonymous(np.dtype(np.float16), "_EqOpaqueHalf") h = gb.Scalar.from_value(1.5, half) @@ -6971,12 +7310,36 @@ def test_udt_numba_error_summary_is_never_empty(): "Something else went wrong" ) assert summarize(ValueError("")) == "Numba could not compile the function for these input types" - # From a lifted op: Numba cannot add two records, so not an array of records either + # From a lifted op: Numba cannot compare two records, so not an array of records either arr_rec = dtypes.register_anonymous( np.dtype(([("sum_p", np.int32), ("sum_q", np.float64)], (2,))), "_SumArrRec" ) - with pytest.raises(UdfParseError, match=r"plus does not work with \(_SumArrRec, _SumArrRec\)"): - binary.plus[arr_rec] + with pytest.raises(UdfParseError, match=r"eq does not work with \(_SumArrRec, _SumArrRec\)"): + binary.eq[arr_rec] + + +@pytest.mark.skipif("not supports_udfs") +def test_udt_records_in_an_array_have_no_arithmetic(): + """A lifted op refuses records in an array with a KeyError that says so. + + It reached Numba, which cannot add two records, and its UdfParseError + quoted Numba's ``add(Record(...), Record(...))``. + """ + rec = [("ria_p", np.int16), ("ria_q", np.float32)] + arr_rec = dtypes.register_anonymous(np.dtype((rec, (2,))), "_RecsInArr") + field = dtypes.register_anonymous(np.dtype([("ria_r", rec, (2,)), ("ria_s", np.float64)])) + match = "a built-in op lifts to the numbers in an array or a record, not to records in" + for udt in [arr_rec, field]: + v = Vector(udt, size=2) + for parent in [binary.plus, binary.min, unary.ainv, unary.abs, monoid.plus]: + with pytest.raises(KeyError, match=match): + parent[udt] + assert udt not in parent + with pytest.raises(KeyError, match=match): + v + v + # Ops that do not touch the values still work. + assert v.ewise_mult(v, binary.first).new().dtype == udt + assert v.apply(unary.identity).new().dtype == udt @pytest.mark.skipif("not supports_udfs") diff --git a/graphblas/tests/test_scalar.py b/graphblas/tests/test_scalar.py index 596de8802..03b416980 100644 --- a/graphblas/tests/test_scalar.py +++ b/graphblas/tests/test_scalar.py @@ -183,6 +183,14 @@ def test_isequal(s): assert s.isequal(t) assert not s.isequal(t, check_dtype=True) assert Scalar.from_value(None, dtype="INT8").isequal(Scalar.from_value(None, dtype="INT16")) + # An int that INT64 cannot hold compares exactly, as with ``==``. + for big in [2**63, 2**64, -(2**63) - 1]: + assert not s.isequal(big) + assert not t.isequal(big) + assert not Scalar(dtypes.INT8).isequal(big) + assert Scalar.from_value(2**63, dtype="UINT64").isequal(2**63) + assert Scalar.from_value(2.0**64).isequal(2**64) + assert not Scalar.from_value(2.0**64).isequal(2**64 - 1) # exactly, as isequal compares @pytest.mark.slow @@ -199,6 +207,12 @@ def test_isclose(): assert not s.isclose(Scalar.from_value(5), check_dtype=True) assert not s.isclose(Scalar.from_value(None, dtype=s.dtype)) assert Scalar.from_value(None, dtype="FP64").isclose(Scalar.from_value(None, dtype="FP32")) + # An int that INT64 cannot hold compares in Python, as in isequal; it raised OverflowError + assert not Scalar.from_value(1, dtype="INT8").isclose(2**63) + assert Scalar.from_value(2.0**64).isclose(2**64 + 1) + assert Scalar.from_value(2**63, dtype="UINT64").isclose(2**63) + assert not Scalar.from_value(1.0).isclose(10**400) + assert not Scalar(dtype="INT8").isclose(2**63) def test_isclose_infinity(): @@ -211,6 +225,10 @@ def test_isclose_infinity(): 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) + # An int that INT64 cannot hold compares in Python, with the same rule + assert not Scalar.from_value(math.inf).isclose(2**70) + assert not Scalar.from_value(-math.inf).isclose(-(2**70)) + assert Scalar.from_value(2.0**70).isclose(2**70) def test_isclose_integers_do_not_wrap(): @@ -224,6 +242,7 @@ def test_isclose_integers_do_not_wrap(): 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) + assert Scalar.from_value(2**64 - 1, dtype="UINT64").isclose(2**64 - 1) def test_isclose_in_the_type_the_dtypes_unify_to(): diff --git a/graphblas/tests/test_ssjit.py b/graphblas/tests/test_ssjit.py index 690e05c6f..8627d80e3 100644 --- a/graphblas/tests/test_ssjit.py +++ b/graphblas/tests/test_ssjit.py @@ -66,26 +66,25 @@ def _jit_works_here(): return _enable_jit_for_udt() -@pytest.fixture(scope="module", autouse=True) -def _setup_jit(): - """Set up the SuiteSparse:GraphBLAS JIT. +def _configure_jit_for_tests(cfg, fix_jit_config): + """Let this module's tests compile JIT kernels where that is cheap; see ``_setup_jit``. Strategy: - 1. _fix_jit_config(): fix conda-baked compiler paths and probe. + 0. Where SuiteSparse compiles through cmake (``jit_use_cmake``, always True + when it is built with MSVC), nothing: as for UDT work + (``_enable_jit_for_udt``), only an ``'on'`` set before the tests, such as + ``GRAPHBLAS_JIT_C_CONTROL=on``, compiles there. Each kernel took about + 10 s through cmake on Windows CI, and this module compiles about 260. + 1. fix_jit_config(): fix conda-baked compiler paths and probe. - Returns True: JIT works, proceed. - Returns False: probe failed. Whatever ``jit_c_control`` SuiteSparse was left at is left untouched. - Returns None: no conda env, try sysconfig instead. 2. Sysconfig fallback: for non-conda installs (pure pip). """ - if _IS_SSGB7: - # SuiteSparse JIT was added in SSGB 8 - yield + if cfg["jit_use_cmake"]: return - - prev = gb.ss.config["jit_c_control"] - - result = _fix_jit_config() + result = fix_jit_config() # ``True``: conda JIT configured and verified. ``False``: probe failed; # don't try sysconfig, since if the conda compiler can't build GraphBLAS # JIT kernels, Python's sysconfig compiler won't either. Both done. @@ -96,12 +95,23 @@ def _setup_jit(): include = sysconfig.get_path("include") libs = sysconfig.get_config_var("LIBS") if cc and cflags and include: - gb.ss.config["jit_c_control"] = "on" - gb.ss.config["jit_c_compiler_name"] = cc - gb.ss.config["jit_c_compiler_flags"] = f"{cflags} -I{include}" + cfg["jit_c_control"] = "on" + cfg["jit_c_compiler_name"] = cc + cfg["jit_c_compiler_flags"] = f"{cflags} -I{include}" if libs: - gb.ss.config["jit_c_libraries"] = libs + cfg["jit_c_libraries"] = libs + +@pytest.fixture(scope="module", autouse=True) +def _setup_jit(): + """Set up the SuiteSparse:GraphBLAS JIT (``_configure_jit_for_tests``) and restore it.""" + if _IS_SSGB7: + # SuiteSparse JIT was added in SSGB 8 + yield + return + + prev = gb.ss.config["jit_c_control"] + _configure_jit_for_tests(gb.ss.config, _fix_jit_config) try: yield finally: @@ -857,6 +867,36 @@ def test_enable_jit_for_udt_re_raises_a_lowered_control(monkeypatch): assert fake["jit_c_control"] == "run" +def test_setup_jit_leaves_the_control_alone_where_suitesparse_uses_cmake(monkeypatch): + """The setup of this module does not set ``'on'`` where SuiteSparse compiles through cmake. + + ``fix_jit_config()`` sets ``'on'`` whenever it finds a compiler to repair + from (on Windows, conda's ``clang.exe`` or the build machine's ``cl.exe``), + so the setup compiled every kernel of this module through cmake there, as + UDT work does not (``_enable_jit_for_udt``). Driven through a stand-in + config, so it runs on every host, whether SuiteSparse was built with MSVC or not. + """ + calls = [] + + def fix_jit_config(): + calls.append(cfg["jit_use_cmake"]) + cfg["jit_c_control"] = "on" + return True + + def no_env(): + calls.append(cfg["jit_use_cmake"]) + + monkeypatch.setattr(sysconfig, "get_config_var", lambda name: f"<{name}>") + for fix in [fix_jit_config, no_env]: + cfg = {"jit_use_cmake": True, "jit_c_control": "run"} + _configure_jit_for_tests(cfg, fix) + assert cfg == {"jit_use_cmake": True, "jit_c_control": "run"} + cfg["jit_use_cmake"] = False + _configure_jit_for_tests(cfg, fix) + assert cfg["jit_c_control"] == "on" + assert calls == [False, False] + + @pytest.mark.skipif("_IS_SSGB7") def test_enable_jit_for_udt_leaves_run_alone_where_suitesparse_uses_cmake(monkeypatch): """UDT work does not raise the default ``'run'`` where SuiteSparse compiles through cmake. @@ -1303,6 +1343,36 @@ def test_jit_udt_name_and_typedef_are_checked(): assert not hasattr(dtypes.ss, "chk_scalar") +@pytest.mark.skipif("_IS_SSGB7") +def test_jit_udt_deserializes_without_dtype(monkeypatch): + """A type from C source is found by its C name, which the blob keeps. + + Such a type has no numpy repr in GrB_NAME, the fallback for a name that does not + resolve, so deserializing without ``dtype=`` raised a bare SyntaxError. + """ + _require_jit_on() + name = "_jit_ser_pair" + np_type = np.dtype([("a", np.int64), ("b", np.float64)]) + dtype = dtypes.ss.register_new( + name, f"typedef struct {{ int64_t a ; double b ; }} {name} ;", np_type=np_type + ) + A = gb.Matrix(dtype, 2, 3) + A[0, 1] = (1, 2.5) + v = Vector(dtype, 3) + v[2] = (3, 4.5) + for obj in [A, v]: + blob = obj.ss.serialize() + got = type(obj).ss.deserialize(blob) + assert got.dtype is dtype + assert got.to_coo()[-1].tolist() == obj.to_coo()[-1].tolist() + # A process that has not registered the type is told to pass it + with monkeypatch.context() as m: + m.delattr(dtypes.ss, name) + with pytest.raises(ValueError, match=f"type name '{name}'; pass it as dtype="): + type(obj).ss.deserialize(blob) + assert type(obj).ss.deserialize(blob, dtype=dtype).dtype is dtype + + def test_jit_unary(v): cdef = "void square (float *z, float *x) { (*z) = (*x) * (*x) ; } ;" if _IS_SSGB7: @@ -2963,6 +3033,25 @@ def test_context_importable(): assert gb.ss.global_context is global_context +@pytest.mark.parametrize("kind", ["binary", "unary"]) +def test_ss_op_module_can_be_the_first_import(kind): + """``graphblas.core.ss.binary`` (or ``unary``) imports as a process's first graphblas import. + + Importing the operator package initializes the operators, which imports + ``graphblas.binary.ss``, which takes ``register_new`` from this module. With + the operator imports at its top, that raised ImportError ("partially + initialized module"), on main too. + """ + script = ( + f"import graphblas.core.ss.{kind} as mod\n" + "import graphblas as gb\n" + "print('package: ' + gb.__file__)\n" + f"print('same: ' + str(gb.{kind}.ss.register_new is mod.register_new))\n" + ) + lines, report = _run_in_fresh_process(script) + assert lines["same"] == "True", report + + def test_fix_from_conda_finds_windows_compilers(tmp_path): """The conda repair looks in ``bin`` and, as Windows conda installs them, ``Library/bin`` with ``.exe`` names; a baked-in Windows path's basename is found too. Nothing there: None. diff --git a/graphblas/tests/test_vector.py b/graphblas/tests/test_vector.py index e7b7e63f7..70b3babd8 100644 --- a/graphblas/tests/test_vector.py +++ b/graphblas/tests/test_vector.py @@ -15,7 +15,6 @@ import graphblas as gb from graphblas import agg, backend, binary, dtypes, indexunary, monoid, select, semiring, unary from graphblas.core import _supports_udfs as supports_udfs -from graphblas.core import lib from graphblas.exceptions import ( DimensionMismatch, DomainMismatch, @@ -2772,12 +2771,9 @@ def test_udt(): vv = Vector.ss.deserialize(v.ss.serialize(), dtype=long_udt) assert v.isequal(vv, check_dtype=True) # The dtype name is too long to embed in the blob; dtype= must be provided. - # Whichever error the caller sees comes from parsing the name that is in - # the blob. Builds with GxB_Serialized_get_String let that SyntaxError - # through; older ones have no second name to try, so they wrap it. - with pytest.raises( - SyntaxError if hasattr(lib, "GxB_Serialized_get_String") else ValueError - ): + # Every build says so (from SuiteSparse 9, the empty GrB_NAME raised a bare + # SyntaxError; before it, there is no second name to try). + with pytest.raises(ValueError, match="pass it as dtype="): Vector.ss.deserialize(v.ss.serialize()) # May be able to look up non-anonymous dtypes by name if their names are too long named_long_dtype = np.dtype([("x", np.bool_), ("y" * 1000, np.float64)], align=False) @@ -3223,6 +3219,28 @@ def test_subarray_dtypes_iso(): np.testing.assert_array_equal(v.to_coo()[1], expected) +# Without SuiteSparse, building a UDT uses dup_op=binary.any, which needs Numba. +@pytest.mark.skipif("not suitesparse and not supports_udfs") +def test_subarray_dtype_values_in_another_array_dtype(): + # Asked for another array dtype of the same shape, the values convert number by + # number. astype to the array dtype repeated each number into an element: shape + # (n, 4, 4), for an iso Vector too. + expected = np.arange(8, dtype=np.int64).reshape(2, 4) + v = Vector.from_coo([0, 1], expected) + iso = Vector("INT64[4]", size=2) + iso[:] = [0, 1, 2, 3] + iso_expected = np.broadcast_to([0, 1, 2, 3], (2, 4)) + for w, values in [(v, expected), (iso, iso_expected)]: + for get in [w.to_dense, lambda w=w, **kw: w.to_coo(**kw)[1]]: + got = get(dtype="FP64[4]") + assert got.dtype == np.float64 + np.testing.assert_array_equal(got, values) + # Leading axes of length 1 do not matter, as in a store + np.testing.assert_array_equal(get(dtype="FP32[1, 4]"), values[:, None, :]) + with pytest.raises(DomainMismatch, match="must have the same shape"): + get(dtype="FP64[2]") + + @pytest.mark.skipif("not supports_udfs") def test_subarray_dtypes(): a = np.arange(3 * 4, dtype=np.int64).reshape(3, 4)