Broadcast return type for MArrays

[ChrisRackauckas]

With the coming changes in v1.0 allowing MArrays to stack-allocate:

https://gist.github.com/ChrisRackauckas/7f3b1664c17ed8350888e217a8721e2f

there isn't a performance reason to prefer SArray over them. Thus, in light of the performance change, I think it may be wise to revisit the choice to make u/t always return an SArray, for both MArrays and SArrays. IMO it seems that simple arithmetic and broadcasting operations should keep mutability since that's a core property of the type, and without it generic algorithms can get messy:

https://stackoverflow.com/questions/46985380/how-do-you-create-an-sarray-marray-with-changed-units

[c42f]

👍 we need to revisit this.

The new changes seem like they've made a huge difference, but it's not all completely straightforward. We should test whether codegen is as good when an MArray is passed to a non-inline function, and various other nontrivial situations.

BTW, I can't foresee MArray ever replacing SArray for non-temporary data structures. Eg Vector{<:SVector} will be inherently more efficient than Vector{<:MVector} because of data layout.

[ChrisRackauckas]

BTW, I can't foresee MArray ever replacing SArray for non-temporary data structures. Eg Vector{<:SVector} will be inherently more efficient than Vector{<:MVector} because of data layout.

I agree, but what I don't think is necessary are the implicit conversions to so greatly favor SArray.

[iyever]

I agree too. For example, when I add two MVectors, I expect the result to have the same type, especially if its an intermediate operation, and I need to mutate it further.

[ChrisRackauckas]

Can this be revisited now that we're squarely in v1.0? We routinely run into this issue because all other mutable array types broadcast to mutable array types, usually retaining the type or choosing a sensible denser container. However, MArrays seem to be incompatible with any arithmetic since they always go to SArrays, so we get PRs to specifically add in type-casting only for MArrays to bring them back to MArrays. I think it's gone too far when downstream packages need to specifically revert the array type changes in order to have compatibility. While I agree users should just use SArrays in most cases, we're trying to set this up for their ease.

[c42f]

It would be great if we could make this work better and I think this would probably be fairly non-breaking (provided performance is still good). I've added this to the list over at #532.

[andyferris]

It should be fine to adjust this here.

Please keep in mind that this is a symptom of a deeper design problem. An AbstractArray can be immutable and AFAICT the only requirements specified by the interface on "functional"-style operations +, map, broadcast, filter, reduce(...; dims = ...) and so-on (as opposed to using similar followed by imperative style map!, broadcast!. filter!, etc) is that an AbstractArray is returned. It needn't be mutable. I'm not even sure if it's documented or not whether that these results can't be lazy and alias the inputs.

So - we can happily adjust this here for the concrete type MArray, but please keep in mind that any code broken by this (that is - ever attempts to mutate the result of a "functional" operation) is, in my opinion, not generic AbstractArray code and liable not to work with arbitrary subtypes of AbstractArray that users might import and use. Of course, that is not a problem when your function signatures dispatch on a concrete type like MArray.

To fix this, we'd need a formal system of array traits (such as mutability - does setindex! work?) and contracts on functions like map that describe which traits are guaranteed to be propagated to the outputs (this is non-trivial - guaranteeing mutability would imply forbidding aliasing/laziness, for example, and I believe both are currently used in the ecosystem). Rather than restricting the output of map, etc, there could new functions to provide a container with a given trait - for example, one that optionally copies an array to a mutable array, but only when necessary. In this case, that function would "upgrade" an SVector to an MVector, but not copy an MVector.

[ChrisRackauckas]

but please keep in mind that any code broken by this (that is - ever attempts to mutate the result of a "functional" operation) is, in my opinion, not generic AbstractArray code and liable not to work with arbitrary subtypes of AbstractArray that users might import and use.

Yes, we know that DiffEq cannot with every AbstractArray, but that's kind of a silly distinction. Technically to be compliant with AbstractArray you have just indexing and cannot assume broadcast, arithmetic, etc. exist, so we are using some subset already. And we put these in a parameterized type, so it requires an AbstractArray whose basic arithmetic is type-stable. Since we are doing math on it and we require efficiency, this doesn't seem to be a big cutoff for our use case. Of all of the input u0 types we have tests for, pretty much all work now because all it really takes is for broadcasted + and * to return the same vector type (note: sparse input doesn't really make sense here). There's only two main cases that are still giving trouble: MArray and Vector{Vector} (though Vector{SVector} is fine). MArray is this issue, and Vector{Vector} and nestings like that are only an issue because there is no recursive broadcast. While having a complete trait system would be amazing, I think the first step in that direction is having common actions act similarly, and right now the issue is that MArray does not work how the vast majority of array types do.

To fix this, we'd need a formal system of array traits (such as mutability - does setindex! work?) and contracts on functions like map that describe which traits are guaranteed to be propagated to the outputs (this is non-trivial - guaranteeing mutability would imply forbidding aliasing/laziness, for example, and I believe both are currently used in the ecosystem). Rather than restricting the output of map, etc, there could new functions to provide a container with a given trait - for example, one that optionally copies an array to a mutable array, but only when necessary. In this case, that function would "upgrade" an SVector to an MVector, but not copy an MVector.

I totally agree that this issue with MArray is just a symptom and not the real issue. But DiffEq is running into issues because we need all of this. I mean, the latest stuff I'm working on ( SciML/OrdinaryDiffEq.jl#549 ) probably looks absolutely bizarre since what is going on has absolutely no documented difference and is as of now just working off of known assumptions carried throughout the package ecosystem (i.e. similar(u,axes(u)) means "the dense thing close to u" in any array type), and it's carrying around "prototype" arrays since we have no way of even describing an abstract array. Here's a set of issues I've opened about parts of generic programming which have undefined AbstractArray handling behavior (which all has a local DiffEq hack):

JuliaLang/julia#21869
JuliaLang/julia#22216
JuliaLang/julia#22218
JuliaLang/julia#22622
JuliaLang/julia#25760
JuliaLang/julia#25107

And that doesn't even cover a lot of the common issues we are now facing, for example the incompatibility of rand! and randn! with AbstractArrays. At least where we are right now is that we are compatible with all of the common array types (well, at least on v0.6. We are still re-building compatibility on v1.0 with a few broadcast cases):

• Array
• Various versions of sparse arrays
• SArray (by having a separate implementation of each algorithm for immutables)
• MArray (we have a workaround for this issue: never use out of place broadcast MArray compatibility SciML/OrdinaryDiffEq.jl#546)
• FieldVector
• ArrayPartition
• VectorOfArray
• DistributedArrays
• GPUArrays (CLArrays and CuArrays)
• AFArrays
• MultiScaleArrays
• LabelledArrays

With this experience we have far more than enough data to know what's required for these generic functions, yet no real documentation of how this is done. I would really enjoy if we could sit down and maybe gather @mbauman to really bang out some prototypes for trait-based array construction, generic array constructors, etc. The broadcast changes were great but were only the first step of what we need.

[ChrisRackauckas]

Finally following through with it. We'll be developing it at https://github.com/JuliaDiffEq/ArrayInterface.jl and upstreaming to Base when we're happy with parts.

[c42f]

Well, we've taken the plunge thanks to #536.