reinterpret SVector as Vector{SVector}

[halleysfifthinc]

Last year, I had code to reinterpret an SVector with length N*3 as a Vector{SVector{3}} like reinterpret(SVector{3, Float64}, SVector{12}(1.0:12.0)). However this now errors:

julia> reinterpret(SVector{3, Float64}, SVector{12}(1.0:12.0))
DimensionMismatch("1:4 is inconsistent with SOneTo{12}")

Stacktrace:
 [1] (::getfield(StaticArrays, Symbol("#errmsg#1")){12})(::UnitRange{Int64}) at ~/.julia/packages/StaticArrays/3KEjZ/src/SOneTo.jl:14
 [2] Type at ~/.julia/packages/StaticArrays/3KEjZ/src/SOneTo.jl:15 [inlined]
 [3] convert at ./range.jl:143 [inlined]
 [4] oftype at ./essentials.jl:370 [inlined]
 [5] axes at ./reinterpretarray.jl:112 [inlined]
 [6] summary at ./show.jl:1863 [inlined]
 [7] show(::IOContext{Base.GenericIOBuffer{Array{UInt8,1}}}, ::MIME{Symbol("text/plain")}, ::Base.ReinterpretArray{SArray{Tuple{3},Float64,1,3},1,Float64,SArray{Tuple{12},Float64,1,12}}) at ./arrayshow.jl:316
 [8] limitstringmime(::MIME{Symbol("text/plain")}, ::Base.ReinterpretArray{SArray{Tuple{3},Float64,1,3},1,Float64,SArray{Tuple{12},Float64,1,12}}) at ~/.julia/packages/IJulia/gI2uA/src/inline.jl:37
 [9] display_mimestring(::MIME{Symbol("text/plain")}, ::Base.ReinterpretArray{SArray{Tuple{3},Float64,1,3},1,Float64,SArray{Tuple{12},Float64,1,12}}) at ~/.julia/packages/IJulia/gI2uA/src/display.jl:67
 [10] display_dict(::Base.ReinterpretArray{SArray{Tuple{3},Float64,1,3},1,Float64,SArray{Tuple{12},Float64,1,12}}) at ~/.julia/packages/IJulia/gI2uA/src/display.jl:96
 [11] #invokelatest#1 at ./essentials.jl:790 [inlined]
 [12] invokelatest at ./essentials.jl:789 [inlined]

The issue being that SOneTo will only convert identical UnitRanges (e.g. SOneTo{4} == 1:4, this is due to a requirement in the constructor at src/SOneTo.jl:12), so the axes are unable to be properly constructed. It seems reasonable to me that if somebody is trying to convert a unit range to an SOneTo, the range of the argument is most likely the intentional one, so relaxing the constructor to something like:

diff --git a/src/SOneTo.jl b/src/SOneTo.jl
index 96c1ad2..3dd29ba 100644
--- a/src/SOneTo.jl
+++ b/src/SOneTo.jl
@@ -9,10 +9,7 @@ end

 SOneTo(n::Int) = SOneTo{n}()
 function SOneTo{n}(r::AbstractUnitRange) where n
-    ((first(r) == 1) & (last(r) == n)) && return SOneTo{n}()
-
-    errmsg(r) = throw(DimensionMismatch("$r is inconsistent with SOneTo{$n}")) # avoid GC frame
-    errmsg(r)
+    (first(r) == 1) && return SOneTo{last(r)}()
 end

 Base.axes(s::SOneTo) = (s,)

Such a change makes my original example work again, and I will note that my example is also quite similar to the example given in #554. I wonder what the original intent behind this restriction was and/or if there are possible negatives to relaxing this? pinging @timholy as the author of that PR

[KristofferC]

Don't you want to do something like

julia> reinterpret(SVector{3, Float64}, [SVector{12}(1.0:12.0),])
4-element reinterpret(SArray{Tuple{3},Float64,1,3}, ::Array{SArray{Tuple{12},Float64,1,12},1}):
 [1.0, 2.0, 3.0]   
 [4.0, 5.0, 6.0]   
 [7.0, 8.0, 9.0]   
 [10.0, 11.0, 12.0]

The current error seems consistent with other basic numerical types:

julia> reinterpret(UInt32, 2.0)
ERROR: bitcast: argument size does not match size of target type

[halleysfifthinc]

Sure, in some sense that could be more valid/consistent. But I disagree for two reasons.

1. My understanding of reinterpret is that for mismatched type sizes (such as the present example, or something like the example you gave for basic numerical types), the second argument needs to be an array (eg reinterpret(UInt32, [2.0]) would work). In my example, the second argument is an (Abstract)array, so wrapping it in an array shouldn't be necessary.

2. Wrapping the second argument in an array would be a copy and/or not be stack allocated anymore, presumably with performance cost as well. In my use case, I am trying to avoid any allocations.

[KristofferC]

argument needs to be an array

It needs to be an Array. The point is that it needs to be on the heap so that it has an identity that the reinterpreted array can reference.

I am trying to avoid any allocations.

You said you wanted to get a Vector{SVector{3}} out which is an allocation not on the stack.

[halleysfifthinc]

It needs to be an Array

If that is the case, then shouldn't that be fixed in the Base code for reinterpret? The only reason my example doesn't currently work is because the axes construction fails. If I fix that, then the example works and there are no allocations. That seems to conflict with your statement that it needs to be on the heap. I've also tested the above patch with incorrectly sized reinterprets, e.g. reinterpret(SVector{5}, SVector{12}(1.0:12.0) which properly fails during the checks in the reinterpret code.

I do want a Vector{SVector{3}} out, but I don't want the allocations; reinterpret allows me to have it both ways. Previously (when this example worked last year), and currently (if I fix the range conversion), I get the intended behavior without any allocations in my function(s).

[andyferris]

AFAICT reinterpret seems only to work on primitive types and as a wraper around Arrays. There used to be a reinterpret of Array element types which used to be more permissive and direct, but this was changed by I think Keno just before Julia 1.0 - nowadays there is a lazy wrapper of the underlying array instead of two arrays which alias each other with different element type (which is a more honest way to describe it LLVM when you enable TBAA optimizations).

Unfortunately you can't reinterpret the bits of either structs or mutable structs directly (I'm not sure if it's possible to do this in the presence of TBAA).

If I understand, what you are asking is to make the ReinterpretArray wrapper work on subtypes of StaticArray? I never realized it was intended for anything other than Array?

[andyferris]

The implementation of ReinterpretArray is fascinating - very low level code dynamically implemented in Julia. It's nice to see GC.@preserve letting you do the C-style pointer/memory stuff safely. It actually makes me wonder if we can specialize the implementation in key cases that I am interested in (Matrix{Float64} <-> Vector{SVector{3,Float64}}).

But yeah, it seems relatively friendly to generic AbstractArray, which is something I didn't know.

@halleysfifthinc I don't think we can relax the constructors of SOneTo to lie, or even for oftype(SOneTo{12}, 1:4) to work, but we might be able to do something. I'm leaning towards a complete implementation of ReinterpretStaticArray.

[halleysfifthinc]

Given that it seems fairly conclusive that the current intent of reinterpret is to be used solely with primitive types, I will reiterate my question: shouldn't the eltype constraints of reinterpret be restricted to primitive types? (Although if primitive types defined outside of base are intended to be supported, it would be more feasible to add that intent/restriction to the docstring?)

+1 for a ReinterpretStaticArray. I had thought of that, but not having the time or expertise, did not pursue it.

For some backstory/clarification, my use case is small numbers (generally < 10) of 3 dimensional vectors/points. When I originally wrote the code, I tested actual Vector{SVector{3}} against reinterpret(SVector{3}, SVector{3*N}(…)) and found better performance with the reinterpret version. I can replicate that even now (with the "fixed" SOneTo conversion) for small N (somewhere between 15 and 20, the Vector{SVector{3}} becomes faster):

julia> @benchmark mean(z) setup=(z = [ SVector{3}(rand(3)) for _ in 1:5 ])
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     4.773 ns (0.00% GC)
  median time:      4.799 ns (0.00% GC)
  mean time:        4.986 ns (0.00% GC)
  maximum time:     25.641 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

julia> @benchmark mean(z) setup=(n = 5; z = reinterpret(SVector{3,Float64}, SVector{n*3}(rand(n*3))))
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     3.618 ns (0.00% GC)
  median time:      3.633 ns (0.00% GC)
  mean time:        3.694 ns (0.00% GC)
  maximum time:     26.129 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

The @code_native output very similar between the two, but I lack the familiarity with assembly to have a better idea as to why the reinterpret version is faster for small N.

Regardless, although I imagine I must have also tested the performance of SVector{N, SVector{3}} last year, (re)testing again shows that this is now clearly the best choice, both in terms of performance and given the consensus that doing a reinterpret shouldn't be possible and/or is a bad ideaᵀᴹ.

julia> @benchmark mean(z) setup=(n=5; z = SVector{n}([ SVector{3}(rand(3)) for _ in 1:n ]))
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     2.273 ns (0.00% GC)
  median time:      2.297 ns (0.00% GC)
  mean time:        2.341 ns (0.00% GC)
  maximum time:     22.108 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

[andyferris]

Given that it seems fairly conclusive that the current intent of reinterpret is to be used solely with primitive types, I will reiterate my question: shouldn't the eltype constraints of reinterpret be restricted to primitive types?

You might have this backwards - I believe that is an implementation detail, not an intention. We'll find out for sure soon enough:

JuliaLang/julia#32660