Reviewing `ismutable` and other mutation / sparsity handling traits

[mcabbott]

Instead of griping on a slack thread & telling people to avoid this package, it was suggested that I make an issue.

I can see that there's some desire to answer the question "can my function work in-place on this array, or should I make a copy?" in a way that allows for Array, StaticArrays, CuArray, SparseArrays, etc. I think that's basically what this package would like to solve, and why it has 1716 dependents by now.

The entire documentation of ismutable appears to be this:

“”"
    ismutable(::Type{T}) -> Bool
Query whether instances of type T are mutable or not, see
https://github.com/JuliaDiffEq/RecursiveArrayTools.jl/issues/19.
“”"

There's also this koan in the readme:

Changed the default of ismutable(array::AbstractArray) = true. We previously defaulted to Base.@pure ismutable(array::AbstractArray) = typeof(array).mutable, but there are a lot of cases where this tends to not work out in a way one would expect.

What one would expect is, however, the whole question. There are quite a few possible answers:

1. true could mean that A[i,j] = x is guaranteed to work (for suitable x, and any i,j in the axes of A). Or perhaps just for some valid indices, not all of them.

2. It could mean that A[i,:] .= x is guaranteed to work (broadcasting not scalar, e.g. CuArrays), or that copyto!(A, B) is guaranteed to work for B of sufficiently similar type (e.g. copyto!(Diagonal(rand(3)), Diagonal(1:3))).

3. However, it could also mean that A[i,j] = x is a desirable way to work. Or that mul!(A, ...) is a good idea. Indexing might work for some sparse arrays, but be very slow. It will also work for self-aliased arrays, but might give weird answers. Maybe you'd like false to warn you off slow paths, or paths where return A isn't obviously the right answer, not just away from errors.

In any of these scenarios, I would guess that the default would be false, to mean "don't try working in-place, fall back to allocating". However:

4. It could be that false means you know there's a good out-of-place way to work. Returned only for SMatrix, etc. Default true will sometimes lead to errors.

5. It could also be that false means you can be sure nobody else will be mutating the array while you work. That's another possible purpose, like a recursive Base.ismutable. In which case the default, on unknown types, has to be true.

Each of these could be useful for someone asking the original "can my function work in-place on this array?" question. But to know whether ismutable answers the question for their particular function, they really need to know which it is.

So I think that the docstring has to very clearly tell us what to expect. What's guaranteed and what's not. What the fallback is and why. Why this isn't Base.ismutable (that one's easy, but still best to mention the name clash). For what purposes you should not use this, and instead use some other function, or just x isa Array.

My guess is that the intended meaning is somewhere between 1 and 2, but I'm not really sure. I presume any array package authors are equally unsure, and hence that how this function gets overloaded will be unreliable. Which is why I recommend against using this package, at present.

Here are some examples where true seems a bit surprising. Are these bugs or are they intended? And have other people spent 5 minutes trying to make up other adversarial examples?

julia> using ArrayInterface, JLArrays, OneHotArrays, ReadOnlyArrays, LinearAlgebra

julia> ArrayInterface.ismutable(UpperTriangular(rand(3,3)))  # can only write some places
true

julia> ArrayInterface.ismutable(Diagonal)  # abstract type, subtypes include Diagonal(1:3)
true

julia> ArrayInterface.ismutable(view(rand(3), [1,1,1]))  # self-aliased, writing will be weird
true

julia> ArrayInterface.ismutable(onehotbatch("abc", 'a':'z'))  # may soon allow setindex!, but a bit like self-aliasing
true

julia> ArrayInterface.ismutable(ReadOnlyArray(rand(3)))  # package designed to block writing
true

julia> ArrayInterface.ismutable(collect("abc")')  # cursed array, can neither read nor write
true

julia> ArrayInterface.can_setindex(jl(rand(3)))  # GPU array, scalar indexing will be bad (same for CuArray)
true

[Tokazama]

I'm not personally aware of what the purpose of ArrayInterface.ismutable at this point. I want to say it's related to some of the desired behavior prior to Base.ismutabletype, but I may be confabulating. Since this package is about arrays (or perhaps in a broader sense collections that have the potential to interact with arrays), I'd use ArrayInterface.can_setindex for most of the desired behavior above.

With regards to which operations are fast, I think it's difficult to have a catch all solution there because you get into the layout of data in addition to the actual ability to mutate individual elements. I'm not sure at what point to consider an operation "fast" when dealing with the added complexity of GPU/CPU, boxed/bits, copy/deepcopy.

[mcabbott]

Well, the only problem is that I don't know what can_setindex does, either. I just thought this gripe-issue was long enough with just one function. So perhaps this should move to a different issue:

help?> ArrayInterface.can_setindex
  │ Warning
  │
  │  The following bindings may be internal; they may change or be removed in future
  │  versions:
  │
  │    •  ArrayInterface.can_setindex

  can_setindex(::Type{T}) -> Bool

  Query whether a type can use setindex!.

julia> x = Diagonal(rand(3));

julia> ArrayInterface.can_setindex(x)
true

julia> for i in eachindex(x)
         x[i] = pi
       end
ERROR: ArgumentError: cannot set off-diagonal entry (2, 1) to a nonzero value (π)

julia> y = ReadOnlyArray(rand(3));

julia> ArrayInterface.can_setindex(y)
true

julia> y[2] = pi
ERROR: CanonicalIndexError: setindex! not defined for ReadOnlyVector{Float64, Vector{Float64}}

julia> z = jl(rand(3));

julia> ArrayInterface.can_setindex(z)
true

julia> z[2] = pi
ERROR: Scalar indexing is disallowed.
Invocation of setindex! resulted in scalar indexing of a GPU array.

[Tokazama]

I'm assuming that the issue for ReadOnlyArray is because we default to true. I'd be in favor of having an opt in approach instead.

As for the Diagonal example, I think we may need to go the same route as isassigned.

[ChrisRackauckas]

No need for inflammatory discussions or titles. Keep it calm and review the code of conduct.

The major interesting cases for that specific set of traits are:

1. Arrays which have the methods defined but shouldn't be used (scalar indexing on GPUs)
2. Arrays that have values which are not all writable (BandedMatrix with A[i,j])
3. Arrays that can freeze/unfreeze
4. Arrays for which you want to avoid A[i,j] at all costs due to performance (SparseMatrixCSC)
5. Arrays that are wrapped in something immutable ( Adjoint{Array} is not Base.mutable)

The original intent of the function was to distinguish for (5) (see SciML/RecursiveArrayTools.jl#19), though it got pretty confusing when later Base.ismutable was added to the language and meant something different 😅, and conflates some of these ideas, so we probably need some kind of breaking change here.

Maybe a more modern way to describe these traits would be to have:

1. is_allowed_setindex(A, I) where I is some index, and this gives true/false for whether index I will error (to avoid having to try/catch)
2. all_allowed_setindex(A) whether one can setindex on all of the indices or just some subset
3. any_allowed_setindex(A) whether any indices can be set
4. fast_scalar_indexing(A) whether scalar indexing is considered to be fast or whether it's considered slow/bad practice for the type (i.e. SparseMatrixCSC or GPUs)

any_allowed_setindex seems almost redundant but there are cases like CuSparseMatrixCSR where it can be expensive to check or find setindexables, so knowing it's allowed somewhere may be a better notion that "ismutable"

I'm not quite sure what to do about sparse matrices though, where setindex is always possible but in many cases that can blow up the non-zeros and so you may want to treat hard zeros as something that's not in the structure.

I can see that there's some desire to answer the question "can my function work in-place on this array, or should I make a copy?" in a way that allows for Array, StaticArrays, CuArray, SparseArrays, etc. I think that's basically what this package would like to solve, and why it has 1716 dependents by now.

This is one small function in the package with like 50. I don't think you can attribute the even a large minority of this package's usage to just this function.

[Tokazama]

If someone want's to tackle a PR for is_allowed_setindex(A, I), all_allowed_setindex, any_allowed_setindex(A), this the PR to Julia where I did this for the same fore isassigned may be a helpful reference JuliaLang/julia#49956. I think reflection on best indexing per array layout opens up a bigger conversation though, since a lot of the work with Chris Elrod on StaticArrayInterface was centered around detecting and propagating this sort of layout info between wrapping layers.