SparseArraysBase.jl release to-do list

[mtfishman]

To-do list for splitting off NDTensors.SparseArraysBase as a separate registered package SparseArraysBase.jl:

• Remove the AnyAbstractSparseArray type union in favor of an @derive macro, similar to Moshi.@derive, the Rust derive attribute for implementing traits, and ArrayLayouts.@layoutmatrix and related macros in ArrayLayouts.jl. This would basically automatically define getindex, map!, etc. as sparse_getindex, sparse_map!, etc. on a specified type or wrapper.
• Reassess which sub-modules of NDTensors.jl that SparseArraysBase depends on, and either remove those dependencies or assess what we need to do to split off those libraries into packages as well, or maybe merge them into SparseArraysBase. For example, right now it relies on:
  • BroadcastMapConversion, which converts broadcast calls to map calls (which is heavily inspired by the broadcasting code logic in Strided.jl). That library is also used in other sub-modules of NDTensors.jl, such as BlockSparseArrays and NamedDimsArrays.
  • TypeParameterAccessors for generically accessing type parameters, in particular it uses functionality for generically getting the type of the parent of a wrapper type. We've been planning to split that off for a while, though I think there are still some type instability issues and interface questions to decide on so I'm not sure how comfortable I am doing that right now.
  • NestedPermutedDimsArrays. This dependency could get turned into a package extension with the sparse array interface getting overloaded with @derive.
• Redesign customization of zero values by storing an array that defines the zero elements at each index, like the design suggested in Proposal: sparse overlay Jutho/SparseArrayKit.jl#13. By default for a sparse array a it could be FillArrays.Zeros(eltype(a), size(a)) (based on FillArrays.Zeros), but for a block sparse array it could be MappedArrays.mappedarray(I -> zeros(eltype(a), blocksizes(a)[I]), CartesianIndices(blocksize(a))) (based on MappedArrays.MappedArray).
• Reassess the minimal interface based on storage_index_to_index, index_to_storage_index, stored_indices, etc. (both the names and functionality). For example, what should be assumed about stored_indices? Should it list the indices in the same order that the corresponding values are stored in the sparse_storage? Should they be expected to have fast (in the DOK case, O(1)) lookup of indices?
  • Assess how the interface of SparseArraysBase corresponds to the current or planned functions in the SparseArrays.jl standard library, such as nnz, nonzerovals (proposed to be changed to storedvals), nonzeroinds (proposed to be changed to storedinds), isstored, etc.
  • Other relevant discussions about sparse iteration APIs: [WIP] generic sparse indices JuliaLang/julia#40103, Sparse matrix interfaces JuliaSparse/SparseArrays.jl#559, add the iteratenz API JuliaSparse/SparseArrays.jl#167, https://github.com/timholy/ArrayIteration.jl, https://github.com/SobhanMP/SparseExtra.jl.
• Introduce a macro @avoid_alloc that changes the behavior of a[B] for block sparse array a so that if the block B doesn't exist in the storage, it returns an unallocated FillArrays.Zeros or UnallocatedArrays.UnallocatedZeros object of the correct size. This can help with optimizing map operations, for example expressions like a[B] += x right now have an unnecessary intermediate allocation, since for the default a::BlockSparseArray (with the default zero constructor), a[B] allocates an array filled with zeros. @avoid_alloc a[B] += x could rewrite that expression as @avoid_alloc BlockSparseArray(blocks(a), z)[B] += x where z = mappedarray(I -> Zeros(eltype(a), blocksizes(a)[I]), CartesianIndices(blocksize(a))), i.e. an array defining the zero elements such that the zero elements of each index are a lazy zero array.
• Create a more robust, general, and efficient system for checking if a function being passed to map/broadcast preserves zero values, and allow customization of that functionality. For example, if the element type is a number/scalar, it could be based on passing a static zero function from VectorInterface.jl, Zeros.jl, StaticNumbers.jl, Static.jl, etc. Alternatively, if the elements are arrays, we could use FillArrays.Zeros. However, these approaches have failure modes, for example if the function being passed isn't generic enough to accept those types, or we can't determine which types or array sizes to use based on the input array.

[lkdvos]

For the dependencies:

• I think the NestedPermutedDimsArray implementation indeed either belongs in a package extension, or maybe it might also be fair to just move it to the NestedPermutedDimsArray module. It's not used in the rest of the SparseArrayBase module, and just consists of an implementation of the interface. I would argue that it would be a bit strange to have SparseArrayBase depend on NestedPermutedDimsArray, rather than the other way around
• For the BroadCastMapConversion, how much different is the implementation from the Strided.jl one? How much hacking does it require to just use the Strided.jl version? From a first glance, it looks like the difference is mostly that Strided only captures StridedView arguments, but it might be possible to either make that more general, or loosen that to AbstractArray since it's already specialized on the broadcast type

[mtfishman]

* I think the NestedPermutedDimsArray implementation indeed either belongs in a package extension, or maybe it might also be fair to just move it to the NestedPermutedDimsArray module. It's not used in the rest of the SparseArrayBase module, and just consists of an implementation of the interface. I would argue that it would be a bit strange to have SparseArrayBase depend on NestedPermutedDimsArray, rather than the other way around

Right, the idea would be to make NestedPermutedDimsArrays.jl into a separate package, and then make a package extension that enables NestedPermutedDimsArray as a wrapper type that is understood by SparseArrayBase.jl. That would tie into the @derive macro that would ideally make it one-liner to register NestedPermutedDimsArray as a wrapper.

* For the `BroadcastMapConversion`, how much different is the implementation from the Strided.jl one? How much hacking does it require to just use the Strided.jl version? From a first glance, it looks like the difference is mostly that Strided only captures StridedView arguments, but it might be possible to either make that more general, or loosen that to `AbstractArray` since it's already specialized on the broadcast type

It's been a while since I wrote it, but I think generalizing from StridedView to AbstractArray is the main substantial difference, besides rewriting it in a way that I found to be easier to understand. I would be ok with that, though I would turn it around and say that this should be functionality that is provided in a separate package since it is more generally useful, so why not have a BroadcastMapConversion.jl package that is used by SparseArraysBase.jl, BlockSparseArrays.jl, Strided.jl, and potentially other packages that are interested in that functionality? I think that one would probably be ready to split off soon, it is pretty small and self-contained (assuming everyone is happy with the interface and implementation, I'd appreciate getting some feedback on it).

EDIT: I think your plan is actually a no-go, since we wouldn't plan to make Strided.jl a dependency of SparseArraysBase.jl and BlockSparseArrays.jl. So splitting it off to be shared seems like the only approach. But maybe what you're saying is that you would prefer to split off the version of the code in Strided.jl rather than the reimplementation in BroadcastMapConversion.jl? As you might imagine, I prefer the BroadcastMapConversion.jl version, since as I said I rewrote it in a way that I thought is easier to understand and has a user-facing interface, but again I'm open to feedback about the code design.

[mtfishman]

Building on the bullet point above about the minimal interface, here is a new interface proposal:

• getstoredindex(a, I...): Get the stored value at index I.
• getunstoredindex(a, I...): Get the unstored/zero value at index I. Defaults to zero(eltype(a)).¹
• eachstoredindex(a): An iterator over the corresponding indices of each stored value (i.e. the same set of indices output by filter(I -> isstored(a, I), eachindex(a))).²
• storedvalues(a): A view of the collection of stored values.
• isstored(a, I...)::Bool: Check if the value at index I is stored or not.
• storedlength(a): The number of stored values. Defaults to length(storedvalues(a)).

I'm not certain that covers everything we need so we'll have to try it out and see. Note that the names are meant to stick as close as possible to the corresponding Base array interface, but using the terminology stored and unstored, i.e. note the similarity to getindex, eachindex, values, isassigned, and length.

Footnotes

1. Note that it could also depend on the trait LinearAlgebra.haszero that would check if the element type has a zero value defined for it just based on the type (i.e. numbers generally do, but non-static arrays generally don't). See also LinearAlgebra.zeroslike for an API for generating zero elements of sparse arrays. ↩

2. See related API proposals: Base.eachstoredindex, ArrayIteration.stored, SparseExtra.iternz. See also the discussions in the SparseArrays.jl standard library here and here. ↩

[lkdvos]

Just checking if I get this right: eachstoredindex(a) gives the indices for the array, so this is equal to filter(I -> isstored(a, I), eachindex(a)), and not equal to eachindex(storedvalues(a)), right?

[mtfishman]

Yes.

[lkdvos]

Okay great. As an aside, it might be useful to have eachstoredindex(a) = eachstoredindex(IndexStyle(a), a), which would allow some more fine-grained control over if these are CartesianIndex entries or simple Ints. I am perfectly okay with defaulting to CartesianIndex, but it's nice to either enforce this, or have a way of selecting one or the other.

[mtfishman]

Sounds like a good idea!