GeezerVison

An Array is a Tensor of stored entries plus a fill value for every other
position, as in pydata/sparse. This first step has the Array class, the
dtypes, creation functions, the element-wise functions, astype and the
data type functions, __array_namespace_info__, basic indexing (ints,
slices, ..., None) and reshape, permute_dims, expand_dims and squeeze.
Element-wise results follow numpy on every position, the fill included;
numpy maps the stored values where no GraphBLAS operator gives its
answer exactly.

…umpy's clip

Every element-wise function is compared with numpy on small arrays of each
dtype, with zero and nonzero fills, NaN, infinities, signed zeros, stored
entries equal to the fill, broadcasting, 0-d arrays and Python scalars.

clip with scalar bounds is numpy's clip, which keeps -0.0 against a bound
of 0.0; with array bounds it is minimum(maximum(x, min), max), as numpy's
clip with array bounds is.

sum, prod, min, max, any and all reduce the stored entries in GraphBLAS
and then account for each fiber's unstored positions from a count of its
entries: k * fill for sum, fill ** k for prod, the fill once for the
others. No count is needed when the fill is the monoid's identity or every
position is stored. min and max propagate NaN, as numpy's do.

mean divides the sum; var and std take two passes (the stored entries'
squared deviations from their fiber's mean in GraphBLAS, plus the fill's
times its count). cumulative_sum and cumulative_prod go through numpy
under the densify limit, since their results are dense. diff subtracts
shifted slices, with prepend and append joined along the axis.

…asic indexing

isdtype takes only this namespace's dtypes (or numpy scalar types), as
numpy's does, not strings.

…ard says

With both fills 0.0, the product at a position only one operand stores is
value * 0.0, which is 0.0 (the fill) only for a finite value with its sign
bit clear: -1.0 * 0.0 is -0.0 and inf * 0.0 is NaN. Each Array caches
whether its stored values are all such, and multiply takes the
intersection (ewise_mult) when both are, the union otherwise.

…-0.0 is present

GraphBLAS's max and min of floats are C's fmax and fmin, which skip NaN
and may return either of 0.0 and -0.0; numpy's maximum and minimum
propagate NaN and order -0.0 below 0.0. Without NaN or -0.0 in either
operand they agree, so an Array now caches a summary of its stored values
(finite, nonnegative, NaN, -0.0), computed once, and maximum and minimum
use GraphBLAS when the summaries and fills allow. min and max reductions
and multiply read the same summary.

Reuse the GraphBLAS Scalars that carry fills and Python scalars (a new
one costs about as much as a small union), skip the broadcast checks for
operands of one shape, and reduce over every axis to a Scalar, whose
count of unstored positions is known without a count reduction.
A 10x10 x + y now costs about 24 us (pydata/sparse: about 110 us; a bare
Tensor union: about 20 us).

An Array that stores every position never shows its fill. Where the
fills of an operation raise (an integer to a negative power) and an
operand is such an Array, the result's fill is never seen either: 0
stands in for it instead of raising an error numpy would not raise.

…ared Matrix once

map_values unloaded each operand's Matrix through a GxB_Container. Two
Tensors that share one Matrix (copy-on-write) made it unload the same
Container twice, and the second read found no values. It now unloads
each Matrix once.

Also drop branches no input reaches: numpy's type resolution already
refuses what they checked.

…every mean

The means of a reduction store exactly the fibers that hold entries, and
every entry lies in one of them, so the deviations are one broadcast
ewise_mult of the entries with the means as they are. Before, the mean of
every fiber was stored densely first (and refused past the densify
limit): var(x, axis=0) of a 1e9 x 1e9 array with a few entries now takes
about a millisecond.

where takes the union of the entries of its three operands (broadcast)
and numpy's where over their aligned values. concat joins the stored
entries in O(entries); the result's fill is the one held by the most
positions, and a part with another fill stores its fill explicitly. stack
is concat of expand_dims.

All cost O(entries of the result): flip slices, moveaxis permutes (O(1)),
and repeat, roll and tile compute the new coordinates of the stored
entries. broadcast_to replicates the entries over the broadcast axes, as
a union does, under the densify limit.

…g_axis

x[key] takes numpy's keys: integer and boolean arrays (as Arrays, numpy
arrays or lists) with ints, slices, ... and None. A boolean index whose
fill is False costs O(its entries). The values are gathered from the
stored entries (Tensor.take or vindex), so a position the array does not
store holds the fill in the result too. numpy's axis order applies:
advanced axes stay in place when adjacent, else go first.

x[key] = value assigns a region by GraphBLAS assign (a basic key, or one
1-d index array among slices, where outer and vectorized indexing agree);
a value whose fill is x's carries its structure, so its unstored
positions delete. Several index arrays scatter: the selected positions
are written, the last of repeated ones winning, as in numpy. Assignment
clears the cached summary of the stored values.

argmax and argmin rank each fiber's stored entries (NaN first, then the
earliest position) and let the fill compete once where the fiber has an
unstored position, at the first such position: numpy's answer in
O(entries log entries). nonzero reads the stored entries when the fill is
zero (every unstored position is nonzero otherwise, so it is dense).
searchsorted searches the stored values of the sorted array and counts
the run of the fill where it sorts.

sort places each fiber's stored values below the fill first, then the run
of the fill (stored values equal to it keep their places in it), then
those above it, NaN last: numpy's stable sort, signed zeros included, in
O(entries log entries). Descending order is flip(sort(flip(x))): values
reversed, equal values in their order, NaN first. argsort's result is a
permutation of every position, so it goes through numpy under the densify
limit.

The unique values are those of the stored entries plus the fill, counted
once per unstored position, so unique_all costs O(entries log entries)
and its inverse_indices keep the input's sparsity (the fill's index is
their fill). NaNs are distinct, as the standard says: a NaN fill with
unstored positions goes through numpy densely, under the limit. isin
checks x1's stored values and fill against x2's.

… extension

Contractions run as GraphBLAS semirings (plus_times, or lor_land for
bool) and treat a zero fill as a structural zero, as scipy.sparse does:
0 * inf at an unstored position does not make NaN. A nonzero finite fill
is the fill times ones plus an array with fill zero, which becomes sums
of the operands broadcast over the result. A fill that is not finite, or
True, goes through numpy densely, under the densify limit.

xp.linalg adds diagonal, trace, vector_norm, the Frobenius, 1 and inf
matrix norms, outer, cross and matrix_power (n >= 0) on the stored
entries, and the decompositions and solvers through numpy densely.
pyproject.toml's codespell list learns 'fro', the standard's name for
the Frobenius norm.

…Array

matmul takes a numpy array operand as an Array, so ndarray @ Array (which
numpy hands to Array.__rmatmul__) works as Array @ Array does.

Matrix.dup holds an n x 1 Matrix by column, as SuiteSparse holds a
Vector, which broke the Tensor's invariant that its Matrix is held by
row (graphblas.array_api's astype of a (2, 1) array failed the check).
The copy is now made into a Matrix created by row.

…products; numpy's clip

By default (graphblas.config array_api_signed_zeros False), a product
with an unstored zero reads 0.0 whatever its sign: multiply of finite
floats with zero fills stores only the intersection, and divide of a
zero-filled numerator by an array with no 0 or NaN stores only the
numerator's entries. A row scaling d[:, None] * A with signed d stored
every position to keep the -0.0s; it now stores A's entries. NaN and inf
stay exact, as does everything else; array_api_signed_zeros True keeps
every signed zero, as numpy.

An intersection of two broadcast operands (an outer product) is now
checked against the densify limit, as a union is.

clip with array bounds is numpy's clip over the union of the entries
(numpy gives the bound on a tie, which minimum(maximum()) did not for
signed zeros).

The tests that compare with numpy bit for bit run in the strict mode.

…nd bool sums

- sum starts from +0.0 as numpy's does, so a sum of -0.0 is 0.0.
- min and max of arrays holding -0.0 (or with a -0.0 fill) reduce with
  numpy's minimum and maximum, which order -0.0 below 0.0; GraphBLAS's
  fmin and fmax take either zero.
- k copies of a complex fill add part by part: inf+0j summed is inf+0j,
  not inf+nanj.
- sum and prod with dtype=bool are logical or and and.
- the counts of unstored positions no longer overflow for fibers of 2**63
  positions or more.
- mean of an empty axis is 0/0 as numpy's: nan, or nan+nanj for complex.
- diff with n=0 ignores prepend and append, as numpy's does.
- asarray of a list of Arrays stacks them (sparse) instead of densifying
  each, and takes numpy arrays of either byte order.

…ll is the identity

With fill 0 (sum) or 1 (product), each fiber holds the fill up to its
first entry, so only the positions from there on are stored. numpy
accumulates one position at a time, so a run of unstored positions adds
the fill once more (0.0 + -0.0 is 0.0): each run becomes one element of
the fiber's sequence, and numpy accumulates the sequences in order, for
numpy's rounding and signed zeros. Other fills still go through numpy
densely.

The diagonal scaling of a broadcast ewise_mult multiplied matrices of
mismatched shapes when the result had a dimension of size 0
(GrB_DIMENSION_MISMATCH, e.g. (2, 1) times (1, 0)). Such a result has no
positions, so it now takes the plain path, as ewise_add does.

…pack of an Array, complex cumsum

- clip with array bounds takes the shape x and its bounds broadcast to,
  as numpy's does (a 0-d x with a bound of shape (0,) is empty).
- from_dlpack of an Array of this namespace needs no copy: copy=False
  returns it, as asarray does.
- cumulative_sum of a complex array with include_initial checked the
  sign of a complex fill with signbit, which takes no complex values.

An Array that stores every position never shows its fill, so the
intersection rules now take its fill as 0: A / d with a dense d (fill 0)
and A * (1 / deg)[:, None] (fill inf) store A's entries instead of every
position. divide by such an array is one broadcast ewise_mult.

Fills are computed with numpy's loop for arrays (1-element arrays), which
numpy's 0-d clip does not share: it broke a tie between 0.0 and -0.0 the
other way.

…r their operands

multiply and divide where one operand's fill is 0 and every value of the
other makes 0 times (or divided by) it a zero now store only the first
operand's entries, whatever the other's fill: A * exp(v)[None, :] with a
sparse v (fill 1) and A / d with d of fill 1. A broadcast second operand
smaller than the first's entries is stored in full first (O(its size))
so that the product is an intersection, not a replication. The strict
mode takes this only where the zeros keep their sign (nonnegative
values), which now includes divide by a positive array.

A fill is computed as numpy computes the dense array: an n-d operand's
fill as a 1-element array and a 0-d operand as 0-d. numpy picks its loop
by those shapes: an array to a 0-d power of 0.5 is sqrt (pow(-inf, 0.5)
is NaN), between two arrays it is pow (inf). The tests seed their
inputs with crc32 of the function's name (hash() of a str varies
between runs), which had made this case show up only in some runs.

matmul, tensordot and vecdot (and their xp.linalg names) take a keyword-only
semiring=, a graphblas Semiring or its name ("min_plus", "max_min", ...), and
the new xp.einsum(subscripts, *operands, semiring=None) runs numpy subscripts
over the Tensor's einsum (plus_times by default, lor_land for bool).

A contraction skips unstored positions, so a semiring other than the default
needs every operand's fill to be the identity of its monoid, and that identity
must annihilate the multiply from both sides (inf for min_plus, -inf for
max_plus, False for lor_land); the result's fill is the identity. The
annihilator is checked once per typed semiring on finite probe values, so
max_times on floats, plus_pair, min_first and integer min_plus raise
ValueError, while unsigned max_times and plus_min pass. The default semiring
keeps the existing rules for nonzero fills; einsum takes numpy's dense path
for them, under the densify limit.

Tensor contractions under a mask now pass an operand stored as the transpose
of the layout the product needs through the transpose descriptor instead of a
transposed copy, so SuiteSparse can use its masked dot product: an SDDMM-shaped
A @ B.T under a mask of 2e5 entries (n = 2e4, k = 64) went from 2.0 s to 1 ms.

|z| of a finite complex z may overflow to inf (|1.5e308 + 1.5e308j|), so abs
of a complex array no longer passes a cached "finite" (or "nan") fact on to
its result: a multiply after it took the intersection and read 0 where numpy
has inf * 0, NaN. Found by review.

…computed once

Four small fixes, found by review:
- The soundness check of a semiring is cached per typed semiring object, not
  by name: two user semirings of one name shared the first one's verdict.
- The SDDMM's dot-product hint divides its saxpy estimate by the number of
  batches (a batched product pairs entries within a batch only). 2000 batches
  of 64 x 64 at density 0.01 under a mask of density 0.5, 8 threads: the
  deferred SDDMM went from 20.4 ms to 7.5 ms (7.8 ms without deferral).
- multiply(P, P) of a pending product computes it once.
- vecdot's docstring says the default semiring, given or not, conjugates x1.

x[0, None] = v with v shaped like x[0, None] raised: the value was sized
without the new axes. For a key of ints, slices, "..." and None, the value's
axes that line up with the new axes (of length 1) are now dropped first. Keys
that mix None with index arrays keep their former behavior. Found by review.

A contraction with an operand whose nonzero fill shows went dense, so
A @ xp.ones(n), einsum("ij,j->i", A, ones) and einsum("ij,->ij", A, 2.0) of a
large sparse A raised DensifyError. An operand with a finite nonzero fill that
is no larger than the largest operand's entries is now stored in full first
(O(its size); an unstored zero of another operand times its finite values is
0 either way), and einsum takes a 0-d operand with a nonzero fill out as a
scalar factor. A non-finite fill still goes dense, as numpy computes it.
Found by review.

An Array could not be pickled: its Tensor holds weak references. Array now
pickles as its stored coordinates and values, shape, dtype and fill, and is
rebuilt by from_coo (explicit fills stay stored; a pending product is
computed first). Found by review.

einsum took a 0-d operand out as a scalar factor of the whole sum; for an
infinite scalar that differs from numpy, which multiplies each term (the
unstored zeros times inf are NaN). A non-finite 0-d operand now takes the
dense path, as before the factoring. The SDDMM's dot hint takes the result's
shape as optional (one batch without it).

…sample

Counting the True values of a fully stored bool result, to decide whether to
drop its False entries, cost as much as the comparison that made it (2e7
values at 1 thread: greater 3.4 ms, with the count 6.4 ms), on every result.
A sample of the leading 2^16 entries now decides: over a third True, the
result is kept as it is (3.3 ms); otherwise one select of the True values
decides, and is the result when at most a quarter of the entries are True.
Dropping still costs the select (2e7 values, 1% True: 22 ms at 1 thread, 4.2 ms
at 8, against 3.2 and 1.4 for the comparison alone). Values never change; a
sample unlike the rest only changes what is stored. Found by review.

…st y's size

Two review findings in the paths that keep products sparse:
- A contraction stored an operand with a finite nonzero fill in full only
  when it was no larger than the other operands' stored entries, so
  A @ xp.ones((n, 8)) of a large sparse A raised DensifyError, and the answer
  could depend on the other operand's explicit zeros (stored in full: no term
  where an unstored zero meets inf; dense: NaN). Such an operand is now stored
  in full whenever its size fits under the densify limit, whatever the others
  store; a non-finite fill still goes dense.
- The one-sided product gathered a broadcast operand's values at every entry
  of the other even when that operand is smaller: A (1e6 x 1e6, 1e7 entries)
  times y[:, None] with y 10% stored took 147 ms and 480 MB more RSS at 8
  threads. It now gathers only when the other has fewer entries than y has
  positions, and otherwise stores y in full for the broadcast product: 6 ms,
  36 MB.

… pattern

einsum("ij,jk,ik->", U, U, U) (triangle counting) built all of U @ U before
the third operand met it. When one operand's indices are exactly what the
other two leave after contracting their shared ones (batch indices
included), the pair is now computed at that operand's stored positions only
(a masked product, with the SDDMM dot hint), times the operand, then summed
into the output. Only for semirings whose multiply distributes over the
monoid (plus_times, lor_land, min_plus, max_plus, min_max, max_min).
Operands holding inf or NaN prune their explicit fills first, as the two-operand path does; a NaN
under a float min or max monoid, a complex operand holding inf, and a
non-finite masking operand under plus_times (inf * (x + y) is not
inf * x + inf * y) take the full einsum. Values do not depend on the path.

…fill selects

Where the condition stores nothing it reads its fill, so the result there is
the operand that fill selects. When that operand needs no broadcast (or is a
scalar, or stores nothing), where now starts from it and writes the
condition's stored positions over it, reading the other operand only at
those positions. SAPS's connected components,
min(where(graph, labels, int_max), axis=1), replicated labels over every
row of graph (n^2 entries; a DensifyError past about 10k nodes); it now
stores graph's entries. Otherwise the union as before.

…ersect

inf + y is inf for every y but -inf and NaN, so when both operands' fills
are inf the result can differ from inf only where both store an entry. Add
(inf or -inf), maximum (inf, or an integer's maximum, or True), minimum (the
opposites) and logical_or (True) now take the intersection when both fills
absorb the op, which never replicates a broadcast operand. An operand that
stores every position counts as holding that fill. Floats need no NaN in
either operand and, for add, no opposite infinity: new cached facts "posinf"
and "neginf" beside "finite", "nan" and "zero". SAPS's Bellman-Ford step,
expand_dims(D, 1) + G with inf fills, wrote D across every column of G.

Such an operand never shows its fill, so matmul, tensordot, vecdot and
einsum with a semiring now relabel its fill as the monoid's identity (O(1), and no
value can tell the fills apart) instead of refusing it, when its dtype holds the identity (an
integer cannot hold inf). A dense distance matrix with no zero, given as an
ndarray, now runs min_plus.

refill(x, fill_value) (an extension) stores the positions that held the old
fill (O(size), under the densify limit; O(1) when x stores every position)
and drops the entries equal to the new fill. A dense distance matrix with
inf for "no edge" and zeros on the diagonal becomes fill inf with its zeros
stored, which a min_plus contraction accepts. The semiring error's hint,
which suggested asarray(x.tensor, fill_value=...) (that turns the positions
holding the old fill into the new one), and asarray's refusal to change an
Array's fill now point at refill.

…s to fill False

equal(deg, 0) of a deg that leaves its few zero degrees unstored has fill
True and stores a False at every other position. When those False entries
are more than half of all positions, the result now takes fill False and
stores only its True positions (its unstored positions and stored True
entries): O(size), under twice its entries. In SAPS's PageRank on a graph
with few dangling nodes, (1/N) * zero_deg then has fill 0 and M stays
sparse (n=2000, 2 dangling: 19,965 entries instead of 3,996,000). A graph
with mostly dangling nodes keeps fill True, which is already the smaller.

matmul, tensordot, vecdot and einsum return an Array that computes on first
use. compute(x) (an extension) computes it at once and returns x, for timing
and for frameworks with a compute step (SAPS's Framework.compute, which had
only the private x._T()). Any other Array is returned as it is.

from_coo raised "Duplicate indices found, must provide dup_op" for a
repeated coordinate; pydata/sparse and scipy.sparse sum such values. It now
sums them (logical or for bool; integers wrap as numpy's sum does). One
value for every coordinate stays stored once (iso) and is built again from
the repeated value only when coordinates repeat.

A 0-d Array stores nothing (its value is its fill), but each one made a new
empty Tensor: a 1 x 1 Matrix plus a per-Matrix format setting, about 6 us,
half of a 0-d multiply. They now take an O(1) copy-on-write copy of one
empty Tensor per dtype (_utils.empty_0d), which a change to any copy does not
reach. A 0-d multiply went from 11.9 to 6.3 us; SAPS's time-stepping loops
make one per operation.

A 0-d operand stores nothing (its value is its fill), so it cannot mask the
pair's product; einsum with one takes the full path. Only an identity-valued
0-d operand reached the masked path before, where the answer was the same.

…other operand

When the operand the condition's fill selects holds only the result's fill
(a scalar, or an array that stores nothing), a condition entry equal to the
fill would store that fill: it is left out, and the other operand is read
only at the remaining entries. where(x > t, x, 0.0) with x of fill 0 stored
every entry of x > t, explicit False ones included (an explicit 0.0 at each
x <= t); it now stores x's entries above t only, as a select would.

…s in GraphBLAS

The previous commit kept only the condition's entries that select the other
operand, but did it in numpy (export the coordinates, gather, rebuild),
single-threaded. When the condition is bool and that operand has the
result's shape, it is now a select of the condition's entries, holding the
operand's fill, overwritten by the operand's entries there (an intersection).
where(x > t, x, 0.0) on 1e7 entries, including the comparison: 0.32 s at 1
thread and 0.058 s at 8 (the numpy path: 0.19 and 0.13 s; before either:
0.33 and 0.068 s, storing twice the entries); the Tensor's select alone takes
0.081 and 0.012 s.

They exported every coordinate and ranked the values with np.lexsort: O(entries log
entries) on one thread, 357x a max along the same axis at 1e7 entries. Now each step is a
GraphBLAS pass over the stored entries, with the fibers as the rows (or columns) of the
stored Matrix: the extreme of each fiber, the entries equal to it (a row or column scaling
by the extremes under eq, with a value at every diagonal position so SuiteSparse runs it as
a scaling), and the smallest position among those (a positional min_firstj or min_secondi
product). The first stored NaN wins over it, as in numpy.

Where the fill beats or ties a fiber's extreme (or is NaN), the fiber's first unstored
position takes part. It is 0 unless the fiber stores position 0, which one positional
product shows; only the fibers that do are read out, to find the first position missing
from their entries. The result stays a Vector until it becomes the output Tensor, and
stores only the positions that are not 0 (its fill).

Where the condition stores nothing its fill selects one operand; when that operand holds
only the fill, the result is the condition's selecting entries holding the other operand. A
broadcast other operand (connected components' where(graph, labels, big)) was read by
exporting the condition's coordinates and gathering in numpy, on one thread. Now the
selecting entries are a select on the condition (truthy ones, or falsy ones under a truthy
fill, for any real dtype), given the other's fill, with its stored values written over by a
broadcast ewise_mult, which the Tensor runs as a diagonal scaling without replicating the
operand. A complex condition still goes through numpy.

A broadcast ewise_mult is a diagonal scaling, but SuiteSparse takes its one-pass scaling
path only when the diagonal has a value at every position; xp.arange(n) stores no 0 (its
fill), and the general product it falls back to is 17x slower at 1e7 entries (0.71 s against
0.041 at 1 thread). The branch is now stored at every position first, when it is no larger
than 4x the condition's selecting entries, and the result is that one product.

…BLAS

where(frontier, k, level), multi-source BFS's level update, has a base (level) that stores
entries. It exported the condition's coordinates, gathered both operands there in numpy,
rebuilt them as a Tensor and took the union with the base: on soc-LiveJournal1 with 8
sources the kernel's memory grew 3.9 GB over its inputs, against 0.85 GB for SAPS's own
python-graphblas framework (which assigns under a mask). Now the condition's selecting
entries are computed as for a base that holds only its fill (a select, then the other
operand's values) and written over the base by one ewise_add: 0.83 GB, and the kernel
takes 2.0 s instead of 3.2. A complex condition still goes through numpy.

The 2-thread cap on a reduce_columnwise into 1e5 or more outputs (603ca79) came from sums:
on FROSTT Uber (4392 x 1.96e6) 4.0 / 4.0 / 17 / 15 ms at 1 / 2 / 4 / 8 threads, and the
same holds now (1e4 x 1e6 with 1e7 entries: 18 / 18 / 11 / 39 ms). MAX and MIN scale
instead, at 4 and 8 threads, with no loss on any shape measured: 5e6 x 5e6 with 1e7 entries
0.122 -> 0.050 s at 8 threads, 1e7 x 1e7 with 1e8 entries 0.66 -> 0.28 s, 1e3 x 1e7
0.100 -> 0.049, Uber 4.1 -> 2.7 ms (JIT off and on alike, SuiteSparse 10.5.1, arm64). So the
cap now applies to the other monoids only.