Referring to a vector in a vector of svectors

New to Julia
1

I am trying to grab a vector within a vector of static vectors. I’ve tried below but it doesn’t seem to work (perhaps since sims is a vector of Svectors and not a vector of vectors?).

sims = Vector{SVector{3,Float64}}(undef, 1000)
data = rand(3,1000)
diff1 = data(1,:) - sims[:][1] 
diff2 = data(2,:) - sims[:][2]

when I do

typeof(sims[:][1])

I get: SVector{3, Float64}

Any ideas?

2

What exactly doesn’t work? What error do you get?

3

You assumed that this would return a vector with 1000 elements that includes all first elements of the SVectors right? That’s not how this works though. You have a vector of SVectors and what you do is retrieve the first element of all elements, which is just the first element, an SVector with 3 elements.

However, because SVectors are packed contiguously in memory you can reinterpret the whole vector as a matrix, I think r = reinterpret(Float64, sims) should create a 3x1000 array. And then you could just do one operation for all your diffs like data .- r. Or slice like r[1, :].

Also you use Matlab slicing notation on data. I think this is where your confusion comes from anyway, because Matlab lets you access struct arrays in a way similar to what you’re trying here.

4

StructArrays provide convenient and performant access to components of all array elements. You can use them instead of regular arrays:

julia> using StructArrays, StaticArrays

julia> sims = StructArray{SVector{3, Float64}}(undef, 10)

# first element, as in a regular array:
julia> sims[1]
3-element SVector{3, Float64} with indices SOneTo(3):
 6.9246147411839e-310
 1.5e-323
 0.0

# first component of each SVector:
julia> sims.:1
10-element Vector{Float64}:
 6.9246147411839e-310
 6.9246148476744e-310
 6.9246154120934e-310
 6.92461485226882e-310
 6.92461491037964e-310
 6.924614910382e-310
 6.9246149103844e-310
 6.9246148522372e-310
 6.92461485224037e-310
 6.92461485224353e-310
5

This does not seem to be the case, but the following is :

r = reshape(reinterpret(Float64, sims), (3,:))
6

Just in case, you could also grab each component using:

getindex.(sims,1)
getindex.(sims,2)
...
3 Likes
7

MATLAB style indexing data(1,:) is (fortunately) not allowed in Julia, use square brackets data[1,:] instead. Also, sims[:][1] doesn’t do what you thought it does. sims[:] is the same as just sims, so this only means sims[1], the first SVector in sims. You cannot subtract 3-element vector from 1000-element vector.

To get first element from all SVectors in sims, use getindex.(sims,1), and to get the second element, use getindex.(sims,2), etc.

This is how you might write that code snippet efficiently in Julia. Notice the use of view and dot broadcasting to avoid allocating extra temporary arrays.

sims = Vector{SVector{3,Float64}}(undef, 1000)
data = rand(3,1000)
diff1 = view(data,1,:) .- getindex.(sims,1) 
diff2 = view(data,2,:) .- getindex.(sims,2)
8

or, more explicitly:

diff1 = [ data[1,i] - sims[i][1] for i in eachindex(sims) ]

ps: you can also write @view(data[1,:]) in the examples above, to avoid copying the data.

9

True, the easiest approach that I use a lot is just tossing a @views macro before the function keyword. This makes the code much cleaner and covers all function scope, I might forget some instances if I write it manually.

@views function ...
end
1 Like
10

As others say sims[:][1] == sims[1], you just need to look at the expression like this: (sims[:])[1]. The expression inside the parentheses is equal to sims itself, except that slicing also makes a copy.

You’re far from the first person I’ve seen to try this approach, so I’m a bit curious about the mental model behind it.

Besides getindex you can use first.(sims) and last.(sims) to get all the first and last elements.

Another thing is that you are accessing sims even though it’s just full of garbage data, initialized with undef. Shouldn’t you put some actual data in there before you start using it?

11

Fyi, the TensorCast package uses similar notation as OP in order to extract the vector components:

using TensorCast
@cast c1[i] := sims[i][1]
@cast c2[i] := sims[i][2]
...
12

Not sure I see the similarity, this is basically just

for i in eachindex(sims) c1[i] = sims[i][1] end

sims[i][1] is exactly what you want, then iterate over i.

But maybe it’s related to array slicing:

c1 = sims[:, 1]

I can see how one might guess that sims[:][1] did something similar. The difference is that [:, 1] is a single action, while [:][1] are two independent operations. I wonder, how do languages where multidimensional indexes are written as [i][j][k] deal with slicing?

13

Thank you everyone. I went the getindex and view route and it worked like a charm.