Add an experimental Array API namespace over Tensor (graphblas.array_api) - #667

…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.

…se product

When an operand of a pending product held inf or NaN, or its products could
overflow, multiply(S, P) checked P's logical size against the densify limit
before computing P whole, so any large sparse product raised DensifyError
where computing without deferral succeeds. It now computes the product (a
sparse product, as without deferral) and multiplies, with no densify check.
Found by review.

…ries

X * y or X / y, where X's fill is 0 and y broadcasts and does not store every
position, stored y in full (O(size of y), within the densify limit) and past
the limit replicated it over its broadcast axes. It now reads y's values at
X's stored coordinates (projected onto y's axes, with y's fill where y
stores nothing): O(entries of X). X of shape (1e4, 1e4, 2) with 4 entries
times exp(W)[:, :, None] (W 1e4 x 1e4 with 1e5 entries): peak RSS 1634 MB to
150 MB, for a result of 4 entries. Found by review.

…nd drop the fill split

Three changes to the contractions, and two small ones, all found by review:
- Explicit fills no longer change values: when an operand holds inf or NaN,
  the operands' stored entries equal to their fill are pruned first, so such a
  value meets an explicit fill as it meets an unstored one, with no term
  ([explicit 0] @ [inf] was NaN, the pruned form 0).
- A stored NaN term propagates through a float min or max monoid, as numpy's
  min and max do; SuiteSparse's fmin and fmax skipped it, so min_plus of
  [[nan, 5]] and [[1], [3]] was 8. The result is computed at once and NaN set
  where a NaN term is: one operand's NaN entries times the others' patterns.
- A nonzero fill no longer splits A into A0 + fill (which turned inf into NaN
  and lost every digit when the fill dwarfs the values): an operand that
  stores every position takes fill 0, any other nonzero fill goes dense.
- The dense path checks the result's size against the densify limit, not
  only the operands'.
- tensordot refuses an axis out of range instead of wrapping it.

|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.