Create a vector of SVectors from a Matrix

General Usage
1

In the MixedModels package one of the structs contains a short, wide matrix for which the primary usage is column-wise. (The struct is VectorFactorReTerm defined in src/modelterms.jl for those who want to play along at home.) To speed up such accesses I create a Vector of Svectors from the columns using

reinterpret(SVector{k,T}, vec(z))

where z is a Matrix{T} of size k by n. This is slightly modified from the suggested form in the documentation for StaticArrays under “Arrays of static arrays”.

Julia v0.6 was okay with this and I could define the field in the struct as Vector{SVector{k,T}}.

In Julia v0.7 I get an error about the type of the field

UBlk: Error During Test at /home/bates/.julia/dev/MixedModels/test/UniformBlockDiagonal.jl:7
  Got exception MethodError(VectorFactorReTerm, (UInt32[0x00000001, 0x00000001, 0x00000001, 0x00000001, 0x00000002, 0x00000002, 0x00000002, 0x00000002, 0x00000003, 0x00000003, 0x00000003, 0x00000003], ["1", "2", "3"], [1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0; -1.0 1.0 -1.0 1.0 -1.0 1.0 -1.0 1.0 -1.0 1.0 -1.0 1.0], [1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0; -1.0 1.0 -1.0 1.0 -1.0 1.0 -1.0 1.0 -1.0 1.0 -1.0 1.0], StaticArrays.SArray{Tuple{2},Float64,1,2}[[1.0, -1.0], [1.0, 1.0], [1.0, -1.0], [1.0, 1.0], [1.0, -1.0], [1.0, 1.0], [1.0, -1.0], [1.0, 1.0], [1.0, -1.0], [1.0, 1.0], [1.0, -1.0], [1.0, 1.0]], :G, ["(Intercept)", "U"], [2], [1.0 0.0; 0.0 1.0], [1, 2, 4]), 0x0000000000006c74) outside of a @test
  MethodError: no method matching VectorFactorReTerm(::Array{UInt32,1}, ::Array{String,1}, ::Array{Float64,2}, ::Array{Float64,2}, ::Base.ReinterpretArray{StaticArrays.SArray{Tuple{2},Float64,1,2},1,Float64,Array{Float64,1}}, ::Symbol, ::Array{String,1}, ::Array{Int64,1}, ::LowerTriangular{Float64,Array{Float64,2}}, ::Array{Int64,1})
  Closest candidates are:
    VectorFactorReTerm(::Array{R,1}, ::Array{String,1}, ::Array{T,2}, ::Array{T,2}, !Matched::Array{StaticArrays.SArray{Tuple{S},T,1,S},1}, ::Symbol, ::Array{String,1}, ::Array{Int64,1}, ::LowerTriangular{T,Array{T,2}}, ::Array{Int64,1}) where {T, R, S} at /home/bates/.julia/dev/MixedModels/src/modelterms.jl:159
    VectorFactorReTerm(!Matched::CategoricalArray{T,1,V,C,U,U1} where U1 where U where C where V where T, !Matched::Array{T,2}, ::Any, ::Any, ::Any) where T at /home/bates/.julia/dev/MixedModels/src/modelterms.jl:171

apparently because the reinterpret is lazy.

So my question is how to instantiate the actual type after the lazy reinterpret. I appreciate that much of the linear algebra manipulations have become lazy and I certainly appreciate that change but I am not yet sure how to say, “I really mean it now - give me one of these instead of a delayed expression formulation.”

2

Have you tried copy( reinterpret(SVector{k,T}, vec(z)) )? This should calculate it immediately!

3

@Seif_Shebl Thanks for the suggestion. I don’t think I want to copy it because I want the same memory locations to be accessible from the original matrix, z, and from the vector of SVectors. Some of the time I treat the structure as a matrix and some of the time I want to treat it as a collection of vectors.

Like most uses of StaticArrays the motivation is that I want some sections of the code to be very fast.

4

Is it feasible to just make your struct parametric?

On Julia 1.1-dev:

julia> x = randn(8,12);

julia> using StaticArrays
[ Info: Precompiling StaticArrays [90137ffa-7385-5640-81b9-e52037218182]

julia> xsv = reinterpret(SVector{8,Float64}, x)
1×12 reinterpret(SArray{Tuple{8},Float64,1,8}, ::Array{Float64,2}):
 [-0.35087, 0.301005, -0.411572, -1.54172, -2.97874, 1.75922, -0.0696284, 0.196542]  …  [-0.42543, -0.951566, 0.273947, 0.653601, -0.954813, 1.81766, -0.671167, -0.00442]

julia> using BenchmarkTools

julia> xsv1 = xsv[1];

julia> @benchmark $xsv[1]
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     4.077 ns (0.00% GC)
  median time:      4.098 ns (0.00% GC)
  mean time:        4.216 ns (0.00% GC)
  maximum time:     15.910 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

julia> @benchmark copy($xsv1)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     2.164 ns (0.00% GC)
  median time:      2.174 ns (0.00% GC)
  mean time:        2.212 ns (0.00% GC)
  maximum time:     7.805 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

Alternatively, if you’re not willing to pay a performance penalty, and can protect the original matrix:

julia> ptr_x = pointer(x);

julia> function getindex_ptr(ptr, ::Type{T}, i) where T
           unsafe_load(Base.unsafe_convert(Ptr{T}, ptr), i)
       end
getindex_ptr (generic function with 1 method)

julia> @benchmark getindex_ptr($ptr_x, SVector{8,Float64}, 1)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     2.404 ns (0.00% GC)
  median time:      2.414 ns (0.00% GC)
  mean time:        2.467 ns (0.00% GC)
  maximum time:     12.073 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

julia> @benchmark getindex_ptr($ptr_x, SVector{8,Float64}, 2)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     2.164 ns (0.00% GC)
  median time:      2.184 ns (0.00% GC)
  mean time:        2.229 ns (0.00% GC)
  maximum time:     19.747 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000
5

You can also consider @Keno’s “trick” of widening the definition of VectorFactorReTerm to accept a reinterpreted matrix. Then your Vector{SVector{k,T}} is concrete and fast and you work with the lazy reinterpreted matrix in MixedModels.

6

Eventually I worked out that what I wanted was collect.

1 Like
7

But collect does the same here as copy?

8

Yes, collect will copy in this case as well. The only way to share memory would be to make the struct parametric.