Convert Vector{SVector} to Matrix

General Usage Performance
1

I have some code where I need to store coordinates x of a mesh. I’m using a vector of points (static vectors), i.e.

using StaticArrays
Point{dim} = SVector{dim, Float64}
x = rand(Point{3}, 10)

but some parts of the code (e.g. using WriteVTK) require the coordinates to be in a matrix instead of a vector. Is there anything wrong with adding the following constructor method to convert the above vector to a matrix

Base.Matrix(v::Vector{SVector{S, T}}) where {S, T} = reduce(hcat, v)

and using

m = Matrix(x)
mT = transpose(Matrix(x))

to get the vector as a matrix and matrix transpose? (EDIT: Yes, because it creates a copy it is not efficient)

Is there a more efficient way to do this? (EDIT: Yes, use reinterpret and reshape as explained below)

Is there any way to do this without making a copy so that modifying elements of matrix m will change the corresponding elements of x? (EDIT: Yes, as above)

2

Checkout the amazing reinterpret:

julia> reinterpret(Float64, x)
30-element reinterpret(Float64, ::Array{SArray{Tuple{3},Float64,1,3},1}):
 0.8552062516276724 
 0.6691039952034261 
...

which you can then reshape. That could maybe be somewhere in the staticarrays doc.

7 Likes
3

Possibly even a convert method? So an array of svectors would be automatically converted to a matrix. But possibly that’d be a bit confusing…

4

Thanks @antoine-levitt, I wasn’t aware of the reinterpret function. This does appear to be much faster for large arrays:

using StaticArrays

Base.Matrix(v::Vector{SVector{S, T}}) where {S, T} = reduce(hcat, v)

eltocols(v::Vector{SVector{dim, T}}) where {dim, T} =
    reshape(reinterpret(Float64, v), dim, :)

eltorows(v::Vector{SVector{dim, T}}) where {dim, T} =
    transpose(reshape(reinterpret(Float64, v), dim, :))

Point{dim} = SVector{dim, Float64}

x = rand(Point{3}, 1000000)

using BenchmarkTools

@benchmark m = Matrix($x)
@benchmark mT = transpose(Matrix($x))
@benchmark m = eltocols($x)
@benchmark mT = eltorows($x)

Results:

julia> @benchmark m = Matrix($x)
BenchmarkTools.Trial: 
  memory estimate:  22.89 MiB
  allocs estimate:  2
  --------------
  minimum time:     10.874 ms (0.00% GC)
  median time:      17.070 ms (0.00% GC)
  mean time:        15.128 ms (4.16% GC)
  maximum time:     53.498 ms (78.78% GC)
  --------------
  samples:          331
  evals/sample:     1

julia> @benchmark mT = transpose(Matrix($x))
BenchmarkTools.Trial: 
  memory estimate:  22.89 MiB
  allocs estimate:  3
  --------------
  minimum time:     10.695 ms (0.00% GC)
  median time:      11.154 ms (0.00% GC)
  mean time:        11.346 ms (2.48% GC)
  maximum time:     52.555 ms (78.71% GC)
  --------------
  samples:          441
  evals/sample:     1

julia> @benchmark m = eltocols($x)
BenchmarkTools.Trial: 
  memory estimate:  64 bytes
  allocs estimate:  2
  --------------
  minimum time:     13.549 ns (0.00% GC)
  median time:      15.117 ns (0.00% GC)
  mean time:        17.357 ns (6.48% GC)
  maximum time:     1.122 μs (97.98% GC)
  --------------
  samples:          10000
  evals/sample:     998

julia> @benchmark mT = eltorows($x)
BenchmarkTools.Trial: 
  memory estimate:  80 bytes
  allocs estimate:  3
  --------------
  minimum time:     16.939 ns (0.00% GC)
  median time:      17.799 ns (0.00% GC)
  mean time:        20.796 ns (8.00% GC)
  maximum time:     1.255 μs (98.28% GC)
  --------------
  samples:          10000
  evals/sample:     998

Moreover, it doesn’t create a copy so the original x vector can be modified by changing the entries of the matrix m.

julia> x = rand(Point{3}, 1000000)
1000000-element Array{SArray{Tuple{3},Float64,1,3},1}:
 [0.4185398323656475, 0.501898669373563, 0.41183968543604865]

julia> m = eltocols(x)
3×1000000 reshape(reinterpret(Float64, ::Array{SArray{Tuple{3},Float64,1,3},1}), 3, 1000000) with eltype Float64:
 0.41854   0.86096   0.762052  0.801235  …  0.293778  0.1493     0.462913
 0.501899  0.964937  0.796707  0.379564     0.932526  0.0662906  0.149578
 0.41184   0.790648  0.52458   0.837944     0.921923  0.477261   0.267126

julia> m[1,1] = 0

julia> x
1000000-element Array{SArray{Tuple{3},Float64,1,3},1}:
 [0.0, 0.501898669373563, 0.41183968543604865]
1 Like
5

It is
https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api.html#Arrays-of-static-arrays-1

4 Likes
6

The doc explains how to convert from a dynamically sized AbstractArray to one of the statically sized array types, and how to reinterpret a dim x N Matrix{T} as a Vector{SVector{dim,T}} but not how to reinterpret a Vector{SVector{dim,T}} as a Matrix{T}.

However, it is still very helpful.

Thank you so much!

7

reinterpret does not give you a Matrix though.

8

Thanks @kristoffer.carlsson, I didn’t think of that. In practice they appear to behave the same (see below) which is what I need. What should I call this “matrix-like” type: 3×1000000 reshape(reinterpret(Float64, ::Array{SArray{Tuple{3},Float64,1,3},1}), 3, 1000000) with eltype Float64? Could it be a an AbstractMatrix? And is there any difference in practice apart from the different type (for dispatch etc.)?

Comparing Matrix(x) with eltocols(x):

julia> m = Matrix(x)
3×1000000 Array{Float64,2}:
 0.468326  0.568942  0.599504  0.569535  …  0.917145   0.440911  0.739385
 0.470099  0.442289  0.811437  0.809661     0.559429   0.42691   0.382616
 0.200765  0.235891  0.985741  0.960975     0.0103481  0.41207   0.507787

julia> r = eltocols(x)
3×1000000 reshape(reinterpret(Float64, ::Array{SArray{Tuple{3},Float64,1,3},1}), 3, 1000000) with eltype Float64:
 0.468326  0.568942  0.599504  0.569535  …  0.917145   0.440911  0.739385
 0.470099  0.442289  0.811437  0.809661     0.559429   0.42691   0.382616
 0.200765  0.235891  0.985741  0.960975     0.0103481  0.41207   0.507787

julia> using LinearAlgebra

julia> norm(m - r)
0.0

julia> norm(m[1,1:10] - r[1,1:10])
0.0
1 Like
9

I guess I just answered my own question:

julia> typeof(m)
Array{Float64,2}

julia> supertype(typeof(m))
DenseArray{Float64,2}

julia> typeof(r)
Base.ReshapedArray{Float64,2,Base.ReinterpretArray{Float64,1,SArray{Tuple{3},Float64,1,3},Array{SArray{Tuple{3},Float64,1,3},1}},Tuple{}}

julia> supertype(typeof(r))
AbstractArray{Float64,2}

julia> typeof(r) <: AbstractMatrix
true
10

Thanks! I was looking for a method to convert a Vector of SVectors into a matrix without allocation, and your technique is the solution.

But as you mentioned, I notice that modifying the resulting ReshapedArray entries modifies the entries of the original SVectors. How is this possible? I thought the entries of SVectors are unmutable. Does this mean that we can always bypass the unmutability of an unmutable type (like SVector in the present example) by reinterpreting it as a mutable type (like ReshapedArray in the present example)?

2 Likes
11

Note you can now also reinterpret(reshape, ...:

julia> x = [@SVector(rand(3)) for _ in 1:5];

julia> reinterpret(reshape, eltype(eltype(x)), x)
3×5 reinterpret(reshape, Float64, ::Vector{SVector{3, Float64}}) with eltype Float64:
 0.12019   0.0829133  0.580475  0.873409  0.710157
 0.97712   0.556687   0.696489  0.419789  0.0425647
 0.506971  0.106399   0.887721  0.116239  0.486331
1 Like
12

They’re elements of a mutable array.
You can always edit the mutable array x in my above example, by replacing one immutable SVector with another immutable SVector.
The reinterpret(reshape array does alias x, so mutating it will also change the contents of the mutable x.

1 Like
13

For example, consider

julia> a = [SVector(1,2,3), SVector(4,5,6)]
2-element Vector{SVector{3, Int64}}:
 [1, 2, 3]
 [4, 5, 6]

julia> A = reinterpret(reshape, Int64, a)
3×2 reinterpret(reshape, Int64, ::Vector{SVector{3, Int64}}) with eltype Int64:
 1  4
 2  5
 3  6

Here, if I set A[1,1] = 0, then a also changes:

julia> A[1,1] = 0
0

julia> a
2-element Vector{SVector{3, Int64}}:
 [0, 2, 3]
 [4, 5, 6]

To me, replacing A[1,1] seems equivalent to replacing the first entry to a[1]. Because a[1] is SVector, being able to replace its first entry seems strange to me.

But you say that replacing A[1,1] is actually equivalent to replacing the entire vector a[1]. I am not sure how replacing one element A[1,1] could be interpreted as replacing the entire a[1] with three entries. Still, I think you are correct, because if we make a an SVector of SVectors, we cannot replace A[1,1]:

julia> a = SVector(SVector(1,2,3), SVector(4,5,6))
2-element SVector{2, SVector{3, Int64}} with indices SOneTo(2):
 [1, 2, 3]
 [4, 5, 6]

julia> A = reinterpret(reshape, Int64, a)
3×2 reinterpret(reshape, Int64, ::SVector{2, SVector{3, Int64}}) with eltype Int64 with indices Base.OneTo(3)×SOneTo(2):
 1  4
 2  5
 3  6

julia> A[1,1] = 0
ERROR: setindex!(::SVector{2, SVector{3, Int64}}, value, ::Int) is not defined.
 Hint: Use `MArray` or `SizedArray` to create a mutable static array
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] setindex!(a::SVector{2, SVector{3, Int64}}, value::SVector{3, Int64}, i::Int64)
   @ StaticArrays ~/Programming/StaticArrays/src/indexing.jl:3
 [3] _setindex_ra!
   @ ./reinterpretarray.jl:648 [inlined]
 [4] setindex!(::Base.ReinterpretArray{Int64, 2, SVector{3, Int64}, SVector{2, SVector{3, Int64}}, true}, ::Int64, ::Int64, ::Int64)
   @ Base ./reinterpretarray.jl:507
 [5] top-level scope
   @ REPL[33]:1
1 Like
14

I think the SVector here is distracting you. We can do this even just with Ints:

julia> b = [1, 3];

julia> B = reinterpret(Bool, b)
16-element reinterpret(Bool, ::Vector{Int64}):
 1
 0
 0
 0
 0
 0
 0
 0
 1
 0
 0
 0
 0
 0
 0
 0

julia> B[1] = false; B[2] = true;

julia> b
2-element Vector{Int64}:
 256
   3

Just like an Int, an SVector{Int} is just a collection of bits. Once we put those bits in an array, the memory of the array is mutable so we can change their values.

1 Like
15

You are not mutating an immutable data structure, if that’s what you are worrying about:

as = SVector(1,2,3)
a = [as, SVector(4,5,6)]
A = reinterpret(reshape, Int64, a)
A[1,1] = 0

Now,

julia> a[1]
3-element SVector{3, Int64} with indices SOneTo(3):
 0
 2
 3

as you know, but as is unchanged.

julia> as
3-element SVector{3, Int64} with indices SOneTo(3):
 1
 2
 3

You have only mutated the mutable variable a.

16

@Mason, this doesn’t explain why we cannot mutate A[1,1] reinterpreted from a that is SVector of SVectors as shown in the last example in my earlier posting above.

17

In your final example, you have an SVector of SVectors as the base memory. That’s all immutable so of course you can’t mutate it, unlike the case where you have an mutable array of SVectors

18

Sure, I’m just just pointing out that my final example contradicts your following claim:

Actually I had the same thought as your above claim, but I was puzzled because my final example didn’t allow me to change the bit array.

19

No it doesnt. Note that I wrote Once we put those bits in an array, the memory of the array is mutable so we can change their values.

SArrays are not mutable. Your example would work if you used an MArray.

20

I feel that the discussion is going a bit tangential, but what did you mean by this?

I thought it meant performing A = reinterpret(reshape, Int64, a). Then it sounded like we could change the contents of A regardless of the type of a, because this operation put a into an array A of bits.