JuliaArrays / ArrayInterface.jl

Designs for new Base array interface primitives, used widely through scientific machine learning (SciML) and other organizations
https://docs.sciml.ai/ArrayInterface/stable/
MIT License
134 stars 37 forks source link

Traits not forwarded for structured matrix wrappers like `Diagonal` #339

Open danielwe opened 2 years ago

danielwe commented 2 years ago

ArrayInterface shows Diagonal{T,<:SVector} as supporting setindex!, which is clearly not true:

using ArrayInterface, ArrayInterfaceStaticArrays, LinearAlgebra, StaticArrays

D = Diagonal(SVector(1.0, 2.0, 3.0))
ArrayInterface.can_setindex(D)
# returns true -- should be false

Should this be fixed somehow? Or is it a misunderstanding of ArrayInterface to think that this would give the correct answer?

danielwe commented 2 years ago

Noticing now that ismutable is forwarded:

julia> ArrayInterface.ismutable(Diagonal([1.0, 2.0, 3.0]))
true

julia> ArrayInterface.ismutable(Diagonal(SVector(1.0, 2.0, 3.0)))
false

Isn't this a bit backwards? Diagonal itself is clearly an immutable type, however the parent array it wraps may or may not be mutable, and setindex! is forwarded. Hence, I might have expected ismutable(::Diagonal) = false (although interpreting mutable as "contains mutable storage" is perhaps more useful), but I definitely expected that can_setindex(::Diagonal) would be forwarded to the parent array---not the other way around.

Tokazama commented 2 years ago

We can't just automatically forward things to parent wrappers because there are cases where the wrapper type changes unknown behavior. Diagonal is also a bit tricky because you can change the value along the diagonal only. We can explicitly fix the example you bring up, but I'm not sure if there's a good solution for all cases.

danielwe commented 2 years ago

Thanks for the response! I understand it's not straightforward, however the current behavior suggests there's been some effort to try to make it work already, albeit incomplete (since forwarding is happening for ismutable, just not can_setindex). Some of the gnarlier cases that come to mind for me are the {Bi,SymTri,Tri}Diagonal types, where the wrapper has multiple parents; however, all these parents at least have the same type.

I haven't spent much time with ArrayInterface, but unless I misunderstand, the idea is that the interface will be implemented by every new array type introduced and supported (ideally this would be done by the developer of the array type, but currently it's mostly you guys). So explicitly fixing Diagonal and the other structured matrix wrappers in Base.LinearAlgebra seems like the correct approach.

Tokazama commented 2 years ago

Right now we end up having to define each method for types in the standard library, so we have to fix it on our end for now. We should probably be more careful about how we forward methods for ismutable, but I'm still not sure what the best approach is for these in-between cases

danielwe commented 2 years ago

Oh yeah, I didn't mean that the fix should live in Base. Let me reword the last sentence in my previous post for clarity:

""" It seems like the correct approach is to implement specific fixes (as opposed to generic automatic forwarding) for Diagonal and the other structured matrix wrappers defined in Base.LinearAlgebra. """

Are you unsure about the desired behavior or just how to implement it without too much hassle/repetition/boilerplate? I don't know your code well enough to have a suggestion about the latter, but the desired behavior seems pretty clear to me: LinearAlgebra's structured matrix wrappers always forward valid setindex! calls to their storage (possibly with a recomputed index), so as far as mutability/setindexability goes, they should inherit traits from the storage type. In particular, I'd argue that as long as you can set some elements (e.g., on the diagonal), can_setindex should return true. (This is reminiscent of GPU arrays without scalar indexing: can_setindex should still return true because it's still allowed to call setindex! with slices. In other words, as long as there exist valid setindex! calls for the type, can_setindex should return true.)

Tokazama commented 2 years ago

That reasoning makes enough sense for those types in LinearAlgebra so we can define it for those. It just can't be generalized to all array types. For example, MappedArrays has some types where the transformation isn't invertible.

Tokazama commented 2 years ago

@danielwe, if you have a clear enough idea of what you think this should do you could make a PR and I'd be happy to review it.