Convert 1D Array to 1D Array of StaticArrays?

General Usage
, ,
1

Hello!

I think this should be fairly simple, suppose I have:

b = collect(1:999)

Where x,y, and z are given as:

x = b[1:3:end]
y = b[2:3:end]
z = b[3:3:end]

How would I then convert / parse this into a Vector of StaticArrays with 333 elements?

Is there a simple way to do it without for loops?

Kind regards

2

Are you looking for something like this?

[b[3*(i-1)+1:3i] for i = 1:333]

You can use StaticArrays if you want.
I would suspect if this for loop can be avoided.

1 Like
3

Yep, that does what I want and is a one liner. I did a small benchmark using BenchmarkTools:

@btime [b[3*(i-1)+1:3i] for i = 1:333]
  16.299 Îźs (669 allocations: 49.64 KiB)

And just for fun reshape

@btime reshape(b,(333,3))
  52.990 ns (2 allocations: 96 bytes)

So I had hoped that it would have been possible to “reshape/convert” b using another way to avoid the allocations. In this specific case it does not matter since both are fast, but for more real scenarios it will start to have an effect on performance for me.

Kind regards

4

You can try using reinterpret if you want to avoid allocations:

reinterpret(SVector{3, eltype(b)}, b)
6 Likes
5

That worked too, giving this performance:

 @btime reinterpret(SVector{3, eltype(b)}, b)
  936.667 ns (10 allocations: 608 bytes)

So a lot better. Do you know though if the type of the reinterpret could become a “headache”? Now it is given as:

Base.ReinterpretArray{SArray{Tuple{3},Int64,1,3},1,Int64,Array{Int64,1}}

Seems like this is the best approach though, other than the return type or am I overthinking?

Kind regards

6

Your benchmark results are skewed due to the fact that b is a non-const global variable (see Manual · BenchmarkTools.jl ). If you fix that, you’ll get more realistic results:

julia> b = collect(1:999);

julia> f_comprehension(b) = [b[3*(i-1)+1:3i] for i = 1:333]
f_comprehension (generic function with 1 method)

julia> f_reshape(b) = reshape(b, (:, 3))
f_reshape (generic function with 1 method)

julia> using StaticArrays

julia> f_reinterpret(b) = reinterpret(SVector{3, eltype(b)}, b)
f_reinterpret (generic function with 1 method)

julia> using BenchmarkTools

julia> @btime f_comprehension($b);
  9.468 Îźs (334 allocations: 39.17 KiB)

julia> @btime f_reshape($b);
  42.198 ns (2 allocations: 96 bytes)

julia> @btime f_reinterpret($b);
  2.964 ns (0 allocations: 0 bytes)

The reinterpret approach is so fast as to be basically free (about 3000 times faster than the comprehension version), and it allocates no memory at all. If a vector of SVectors is what you want, then this is definitely the way to go.

As for the return type, it should be fine. You’ll need to make sure that the functions you pass that to expect an AbstractVector or an AbstractArray, but that’s good practice anyway.

5 Likes
7

I see thanks, marked reinterpret as the answer :slight_smile:

And thanks for the hints regarding using AbstractArrays etc., I did not think of that usually I would define numeric type too.

Kind regards

8

It is just ugly :slight_smile: But apparently there is no “solution” to the aesthetic problem.

1 Like
9

I would like to add, because I just faced this, that for the reinterpretation of a 2D array to a vector of static vectors to actually work in function accepting AbstractVectors, one need to add a vec to the conversion:

julia> f(x::AbstractVector{<:AbstractVector}) = "this is a vector of vectors"
f (generic function with 1 method)

julia> x = rand(3,5)
3×5 Matrix{Float64}:
 0.33453   0.209665  0.0140571  0.385981  0.0500048
 0.922044  0.549503  0.826726   0.875412  0.200793
 0.418103  0.830735  0.913105   0.619941  0.740944

julia> y = vec(reinterpret(SVector{3,Float64},x)) # works in f
5-element reshape(reinterpret(SVector{3, Float64}, ::Matrix{Float64}), 5) with eltype SVector{3, Float64}:
 [0.3345302406715238, 0.9220439220026388, 0.4181031420691954]
 [0.2096650742720747, 0.5495028980965941, 0.8307349539720923]
 [0.014057112572854358, 0.8267260218514032, 0.9131049059956979]
 [0.3859812308760291, 0.8754122419901476, 0.6199408169713525]
 [0.050004821335951855, 0.2007929870776446, 0.7409443625689442]

julia> f(y)
"this is a vector of vectors"

julia> z = reinterpret(SVector{3,Float64},x) # does not work in f
1×5 reinterpret(SVector{3, Float64}, ::Matrix{Float64}):
 [0.33453, 0.922044, 0.418103]  [0.209665, 0.549503, 0.830735]  …  [0.385981, 0.875412, 0.619941]  [0.0500048, 0.200793, 0.740944]

julia> f(z)
ERROR: MethodError: no method matching f(::Base.ReinterpretArray{SVector{3, Float64}, 2, Float64, Matrix{Float64}, false})
Closest candidates are:
  f(::AbstractVector{var"#s2"} where var"#s2"<:(AbstractVector{T} where T)) at REPL[24]:1
Stacktrace:
 [1] top-level scope
   @ REPL[29]:1

edit: but noticed now that that is not needed for the 1D array which is what the OP asked, because the array was already a vector and the reinterpreted array carries the type hierarchy of the original array (is that correct to say?)

1 Like
10

What is the purpose of the ‘ugly’ type signatures for reshape and ‘reinterpret’? Is there any difference in layout between

and Array{SVector{3, Int64},1}?

I thought that the only thing that reshape and reinterpret did was to change the type tag on the underlying data. If so, why not simplify the type?

I seem to vaguely remember that is used to be different.

1 Like
11

Is that better than:

x = rand(3,5)
y = SVector{5}(reinterpret(SVector{3,Float64},x))

?

12

In my MWE I didn’t mean to imply that 5 was typical (small) dimension. Usually I’m talking about particle simulations, and this dimension is on the tenths-of-thousands. So making that a static array is not possible.

1 Like