Rewrite SLVector

[MSeeker1340]

This is based on @YingboMa's new SLVector implementation and the focus is on fixing the general conversion behavior of an SLVector: convert(AbstractVector{S}, x) and the broadcasted S.(x)/convert.(S, x).

My attempt is to define a UnionAll type SLVector and have special methods defined on that. For example, I overwrote AbstractVector{S}(x::SLVector) (which is called by convert(AbstractVector{S}, x) by default) because the default implementation uses copyto! on a similar vector, but since SVectors are immutable this is problematic, so I completely overhauled the behavior.

(On a side note, the fact the AbstractVector{S}(x) uses copyto! has prompted StaticArrays to implement similar(x::SVector) to return an MVector (so it can be copied to). I don't really like this but I understand that they're trying to keep the behavior consistent with Base. Similar story goes to broadcast handling, which I also overwrote.)

The broadcasting implementation is broken at the moment: for some reason calling broadcast as a function yields the correct result, but using the dot syntax doesn't.

julia> ABC = @SLVector Int (:a, :b, :c)
#3 (generic function with 1 method)

julia> x = ABC(1,2,3);

julia> y = ABC(4,5,6);

julia> x .+ y
3-element LVector{Int64,Array{Int64,1},(:a, :b, :c)}:
 5
 7
 9

julia> broadcast(+, x, y)
3-element LVector{Int64,SArray{Tuple{3},Int64,1,3},(:a, :b, :c)}:
 5
 7
 9

julia> Float64.(x)
3-element LVector{Float64,Array{Float64,1},(:a, :b, :c)}:
 1.0
 2.0
 3.0

julia> broadcast(Float64, x)
3-element LVector{Float64,SArray{Tuple{3},Float64,1,3},(:a, :b, :c)}:
 1.0
 2.0
 3.0

This seems really bizzare to me as the Julia docs state that f.(x) and broadcast(f, x) are equivalent (https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting-1).

[MSeeker1340]

There's also an error in the test script diffeq.jl when using SLVector in an ODEProblem. The reason is that in initialization, integrator.k is now correctly an SLVector (previously because of conversion problems it's just a regular LVector) and hence immutable, but we're trying to set its fields (https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/perform_step/low_order_rk_perform_step.jl#L528)

We can either go through OrdinaryDiffEq and make sure all out-of-place do not modify the vectors in any way (which could be a big endeavor) or just use MVector instead of SVector. In light of this (JuliaArrays/StaticArrays.jl#327) and @YingboMa's latest PR to StaticArrays I think MVector should be the way to go.

[ChrisRackauckas]

(On a side note, the fact the AbstractVector{S}(x) uses copyto! has prompted StaticArrays to implement similar(x::SVector) to return an MVector (so it can be copied to). I don't really like this but I understand that they're trying to keep the behavior consistent with Base. Similar story goes to broadcast handling, which I also overwrote.)

Similar is defined as returning a mutable type. See JuliaLang/julia#22218 . This is one of the issues referenced at https://github.com/JuliaDiffEq/ArrayInterface.jl .

[ChrisRackauckas]

This seems really bizzare to me as the Julia docs state that f.(x) and broadcast(f, x) are equivalent (https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting-1).

That's wrong. f.(x) is equivalent to broadcast((_x)->f(_x),x), so if someone dispatched on broadcast(::typeof(f),x) then you'll get a different method.

[ChrisRackauckas]

There's also an error in the test script diffeq.jl when using SLVector in an ODEProblem. The reason is that in initialization, integrator.k is now correctly an SLVector (previously because of conversion problems it's just a regular LVector) and hence immutable, but we're trying to set its fields (https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/perform_step/low_order_rk_perform_step.jl#L528)

That's why it integrator.k was mutable. I agree that it's better to just solve this problem but that's a larger issue to handle.

[MSeeker1340]

That's wrong. f.(x) is equivalent to broadcast((_x)->f(_x),x), so if someone dispatched on broadcast(::typeof(f),x) then you'll get a different method.

The behavior is different even when f is a generic user-defined function:

julia> function func(x)           
           println("here")
           x
       end
func (generic function with 1 method)

julia> func.(x)
here
here
here
3-element LVector{Int64,Array{Int64,1},(:a, :b, :c)}:
 1
 2
 3

julia> broadcast(func, x)
here
here
here
3-element LVector{Int64,SArray{Tuple{3},Int64,1,3},(:a, :b, :c)}:
 1
 2
 3

Could the issue be inlining then? Anyways I have come to understand that overload broadcast directly is probably not a good idea (e.g. the error in the test script when I try to do mixed type addition broadcasting), however the broadcast handling interface recommended by the language (broadcast style and such) is also hard to adapt to SLVector since the broadcasting machinery, like AbstractVector{T}, defaults to allocate cache then copy, which is not possible if the array is mutable.

I think this is another place where switching to MArray makes the most sense, as that would make things significantly easier (I only need to modify the exisiting broadcasting handling functions by Yingbo so that they use similar(x.__x, S) as the output cache instead of hard-coded to use Vector{S}.

[YingboMa]

julia> Meta.@lower x .+ x
:($(Expr(:thunk, CodeInfo(
 1 ─ %1 = (Base.broadcasted)(+, x, x)                                                                                                                       │
 │   %2 = (Base.materialize)(%1)                                                                                                                            │
 └──      return %2                                                                                                                                         │
))))

[YingboMa]

SArray doesn't use the similar(bc::Broadcast.Broadcasted{...}, ::Type{ElType}) machinery.
https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/broadcast.jl#L93

That is the reason that they are different.

[ChrisRackauckas]

No, that's just an internal function from the broadcast handling:

https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/broadcast.jl#L6-L38

@andyferris do you know why that would happen?

[YingboMa]

https://github.com/JuliaLang/julia/blob/master/base/broadcast.jl#L767

Normally, copy calls similar(bc::Broadcast.Broadcasted{...}, ::Type{ElType}) and performs copyto!, but SArray is immutable and it calls _broadcast. That is what I mean by it doesn't use the similar machinery.

https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/broadcast.jl#L24

[YingboMa]

Maybe we can overload copy! and copyto! for SLVectors to avoid this problem.

[MSeeker1340]

Oops... wasn't looking carefully when doing the diffeq.jl test error analysis.

There's also an error in the test script diffeq.jl when using SLVector in an ODEProblem. The reason is that in initialization, integrator.k is now correctly an SLVector (previously because of conversion problems it's just a regular LVector) and hence immutable, but we're trying to set its fields (https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/perform_step/low_order_rk_perform_step.jl#L528)

What is actually happening here is: integrator.k is a regular array of SLVector (as it should be). The array assignment itself is not the problem: it's the use of zero, which defaults (again) to creating a similar copy and filling it zeros, thus the setindex! error in the stack trace.

SVector has zero overloaded so using it is fine, we only need to do the same for SLVector here.

But this reminded me that we should go through StaticArrays and copy all the overloaded methods also to here. Julia base is really not friendly towards immutable arrays :( kinda wish it includes a immutable_array trait that the base functions can dispatch on instead of treating all AbstractArray as mutable, or just introduce an AbstractImmutableArray type.

[ChrisRackauckas]

But this reminded me that we should go through StaticArrays and copy all the overloaded methods also to here. Julia base is really not friendly towards immutable arrays :( kinda wish it includes a immutable_array trait that the base functions can dispatch on instead of treating all AbstractArray as mutable, or just introduce an AbstractImmutableArray type.

See and use ArrayInterface.jl

SVector has zero overloaded so using it is fine, we only need to do the same for SLVector here.

For now we can Lazy.@forward a lot of stuff. Though I wonder if having SLVector be a different struct which is <:StaticArray could be better/easier.

[MSeeker1340]

Though I wonder if having SLVector be a different struct which is <:StaticArray could be better/easier.

I can see the appeal to this. On the other hand I really like Yingbo's proposed solution since it treats LVector and SLVector uniformly... but I guess practicability comes first.

For now we can Lazy.@forward a lot of stuff.

What's this Lazy.@forward?

[MSeeker1340]

Closed in favor of #19.