diff --git a/docs/user_guide/collections.rst b/docs/user_guide/collections.rst index de7469c6d..85af39125 100644 --- a/docs/user_guide/collections.rst +++ b/docs/user_guide/collections.rst @@ -116,6 +116,32 @@ SuiteSparse:GraphBLAS has very good defaults which work for most use cases witho When manual tuning is required, look in the :ref:`api_reference` for details about the ``.ss`` namespace. +One pattern worth tuning is a Vector written by a masked assign and read back as a mask in every round, +such as the levels of a breadth-first search: + +.. code-block:: python + + levels = gb.Vector(int, A.nrows) + levels[source] = 0 + levels.ss.config["sparsity_control"] = "bitmap" + frontier = gb.Vector(bool, A.nrows) + frontier[source] = True + depth = 1 + while True: + frontier(~levels.S, replace=True) << frontier.vxm(A, gb.semiring.any_pair) + if frontier.nvals == 0: + break + levels(frontier.S) << depth + depth += 1 + +While ``levels`` is sparse, the masked assign leaves its new entries pending, and reading ``~levels.S`` +as the next round's mask makes SuiteSparse:GraphBLAS merge them into the Vector first. A bitmap Vector +takes the assign in place, so no round waits. The cost is memory once the Vector fills: a bitmap keeps +a byte per position beside the values, and it cannot become the smaller full format. + +Set the option once the Vector has an entry, as above. A masked assign into an empty Vector replaces +its ``sparsity_control`` with the mask's (SuiteSparse:GraphBLAS 10.5). + Collection Methods ------------------ diff --git a/graphblas/core/matrix.py b/graphblas/core/matrix.py index 508a7988a..717a16ead 100644 --- a/graphblas/core/matrix.py +++ b/graphblas/core/matrix.py @@ -50,6 +50,7 @@ from .utils import ( _CArray, _Pointer, + _repeated_value, check_element_shape, class_property, get_order, @@ -1035,6 +1036,8 @@ def from_coo( Column indices. values : list or np.ndarray or scalar, default 1.0 List of values. If a scalar is provided, all values will be set to this single value. + A 1-d array that repeats one value with stride 0, as ``np.broadcast_to(x, n)`` + makes, is stored once (iso) as a scalar is, unless ``dup_op`` is given. dtype : Data type of the Matrix. If not provided, the values will be inspected to choose an appropriate dtype. @@ -1063,6 +1066,14 @@ def from_coo( """ rows = ints_to_numpy_buffer(rows, np.uint64, name="row indices") columns = ints_to_numpy_buffer(columns, np.uint64, name="column indices") + # n copies of one value (as from np.broadcast_to) build iso. Unlike a + # scalar value, they must still number nvals, and duplicates still raise. + repeated = None + if dup_op is None and backend == "suitesparse": + repeated = _repeated_value(values, dtype) + if repeated is not None: + n = values.size + values = repeated values, dtype = values_to_numpy_buffer(values, dtype, subarray_after=1) # Compute nrows and ncols if not provided if nrows is None: @@ -1075,7 +1086,16 @@ def from_coo( ncols = int(columns.max()) + 1 # Create the new matrix C = cls(dtype, nrows, ncols, name=name) - if values.ndim == 0: + if repeated is not None: + if rows.size != n or columns.size != n: + raise ValueError( + "`rows` and `columns` and `values` lengths must match: " + f"{rows.size}, {columns.size}, {n}" + ) + C.ss.build_scalar(rows, columns, values.tolist()) + if C._nvals < n: + raise ValueError("Duplicate indices found, must provide `dup_op` BinaryOp") + elif values.ndim == 0: if dup_op is not None: raise ValueError( "dup_op must be None if values is a scalar so that all " diff --git a/graphblas/core/utils.py b/graphblas/core/utils.py index ddad1d327..86105f960 100644 --- a/graphblas/core/utils.py +++ b/graphblas/core/utils.py @@ -171,6 +171,24 @@ def _distinct_buffers(*arrays): return rv +def _repeated_value(values, dtype): + """Return the one value a stride-0 1-d array repeats, as a 0-d array, or None. + + ``np.broadcast_to(x, n)`` spells n copies of ``x`` without making them. + Values of a UDT (from ``values`` or ``dtype``) give None. + """ + if ( + type(values) is np.ndarray + and values.ndim == 1 + and values.size > 0 + and values.strides == (0,) + and values.dtype.kind in "biufc" + and (dtype is None or not lookup_dtype(dtype)._is_udt) + ): + return values[0, ...] + return None + + def _get_subdtype(dtype): while dtype.subdtype is not None: dtype = dtype.subdtype[0] diff --git a/graphblas/core/vector.py b/graphblas/core/vector.py index 49325b34e..451be6ccb 100644 --- a/graphblas/core/vector.py +++ b/graphblas/core/vector.py @@ -40,6 +40,7 @@ from .utils import ( _CArray, _Pointer, + _repeated_value, check_element_shape, class_property, ints_to_numpy_buffer, @@ -825,6 +826,8 @@ def from_coo(cls, indices, values=1.0, dtype=None, *, size=None, dup_op=None, na Vector indices. values : list or np.ndarray or scalar, default 1.0 List of values. If a scalar is provided, all values will be set to this single value. + A 1-d array that repeats one value with stride 0, as ``np.broadcast_to(x, n)`` + makes, is stored once (iso) as a scalar is, unless ``dup_op`` is given. dtype : Data type of the Vector. If not provided, the values will be inspected to choose an appropriate dtype. @@ -850,6 +853,14 @@ def from_coo(cls, indices, values=1.0, dtype=None, *, size=None, dup_op=None, na """ indices = ints_to_numpy_buffer(indices, np.uint64, name="indices") + # n copies of one value (as from np.broadcast_to) build iso. Unlike a + # scalar value, they must still number nvals, and duplicates still raise. + repeated = None + if dup_op is None and backend == "suitesparse": + repeated = _repeated_value(values, dtype) + if repeated is not None: + n = values.size + values = repeated values, dtype = values_to_numpy_buffer(values, dtype, subarray_after=1) # Compute size if not provided if size is None: @@ -858,7 +869,15 @@ def from_coo(cls, indices, values=1.0, dtype=None, *, size=None, dup_op=None, na size = int(indices.max()) + 1 # Create the new vector w = cls(dtype, size, name=name) - if values.ndim == 0: + if repeated is not None: + if indices.size != n: + raise ValueError( + f"`indices` and `values` lengths must match: {indices.size} != {n}" + ) + w.ss.build_scalar(indices, values.tolist()) + if w._nvals < n: + raise ValueError("Duplicate indices found, must provide `dup_op` BinaryOp") + elif values.ndim == 0: if dup_op is not None: raise ValueError( "dup_op must be None if values is a scalar so that all " diff --git a/graphblas/tests/test_matrix.py b/graphblas/tests/test_matrix.py index c1efcdb70..e6a253de6 100644 --- a/graphblas/tests/test_matrix.py +++ b/graphblas/tests/test_matrix.py @@ -179,6 +179,33 @@ def test_from_coo_scalar(): C.ss.iso_value +@pytest.mark.parametrize("value", [np.float64(1.5), np.int32(-3), [True], np.uint64(2**64 - 1)]) +@pytest.mark.parametrize("dtype", [None, "FP32"]) +def test_from_coo_broadcast_values(value, dtype): + # n copies of one value that were never made (stride 0) build iso, and + # otherwise give what the same values as an ordinary array give. The dtype + # comes from the array: np.int32 values stay INT32, where np.int32(-3) as + # the values gives INT64. + values = np.broadcast_to(value, 3) + C = Matrix.from_coo([0, 1, 3], [1, 1, 2], values, dtype) + expected = Matrix.from_coo([0, 1, 3], [1, 1, 2], values.copy(), dtype) + assert C.isequal(expected, check_dtype=True) + if suitesparse: + assert C.ss.is_iso + + values = np.broadcast_to(np.float64(1), 3) + with pytest.raises(ValueError, match="Duplicate indices found"): + Matrix.from_coo([0, 1, 0], [1, 1, 1], values) + with pytest.raises(ValueError, match="lengths must match: 2, 2, 3"): + Matrix.from_coo([0, 1], [1, 1], values) + # dup_op combines duplicates as it does any values, so the result isn't iso + C = Matrix.from_coo([0, 1, 0], [1, 1, 1], values, dup_op=binary.plus) + assert C.isequal(Matrix.from_coo([0, 1], [1, 1], [2.0, 1.0]), check_dtype=True) + C = Matrix.from_coo([], [], np.broadcast_to(np.float64(1), 0), nrows=2, ncols=2) + assert C.nvals == 0 + assert C.dtype == dtypes.FP64 + + @pytest.mark.skipif("not suitesparse") def test_ss_iso_value_of_an_empty_iso_matrix(): # SuiteSparse marks a Matrix built from no tuples iso, and may mark other empty diff --git a/graphblas/tests/test_vector.py b/graphblas/tests/test_vector.py index 0bfa80b21..69849c47b 100644 --- a/graphblas/tests/test_vector.py +++ b/graphblas/tests/test_vector.py @@ -249,6 +249,27 @@ def test_update_with_a_literal_raises_graphblas_errors(): v(Vector(bool, 3).S) << 1 +@pytest.mark.parametrize("value", [np.float64(1.5), np.int32(-3), [True], np.uint64(2**64 - 1)]) +@pytest.mark.parametrize("dtype", [None, "FP32"]) +def test_from_coo_broadcast_values(value, dtype): + # As for Matrix.from_coo: stride-0 values build iso, and otherwise give + # what the same values as an ordinary array give + values = np.broadcast_to(value, 3) + v = Vector.from_coo([0, 1, 3], values, dtype) + expected = Vector.from_coo([0, 1, 3], values.copy(), dtype) + assert v.isequal(expected, check_dtype=True) + if suitesparse: + assert v.ss.is_iso + + values = np.broadcast_to(np.float64(1), 3) + with pytest.raises(ValueError, match="Duplicate indices found"): + Vector.from_coo([0, 1, 0], values) + with pytest.raises(ValueError, match="lengths must match: 2 != 3"): + Vector.from_coo([0, 1], values) + v = Vector.from_coo([0, 1, 0], values, dup_op=binary.plus) + assert v.isequal(Vector.from_coo([0, 1], [2.0, 1.0]), check_dtype=True) + + def test_from_coo_default_values(): u = Vector.from_coo([1, 3]) v = Vector.from_coo([1, 3], [1.0, 1.0])