Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add a note about reinterpret's memory layout #199

Merged
merged 8 commits into from
Aug 27, 2021
Merged
Changes from 5 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 67 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -342,3 +342,70 @@ julia> s
Foo(44, "d")
Foo(55, "e")
```

## Advanced: StructArrays versus struct-of-arrays layout in higher-dimensional array

Regular arrays of structs can sometimes be reinterpreted as arrays of primitive values with an added
initial dimension.

```julia
julia> v = [1.0+3im, 2.0-im]
2-element Vector{ComplexF64}:
1.0 + 3.0im
2.0 - 1.0im

julia> reinterpret(reshape, Float64, v)
2×2 reinterpret(reshape, Float64, ::Vector{ComplexF64}) with eltype Float64:
1.0 2.0
3.0 -1.0
```

However, the situation is more complex for the `StructArray` format, where `s = StructArray(v)` is
stored as two separate `Vector{Float64}`. `reinterpret` on `StructArray` returns an
"array-of-structs" layout, as the reinterpretation works element-wise:

```julia
julia> s = StructArray([1.0+3im, 2.0-im])
2-element StructArray(::Vector{Float64}, ::Vector{Float64}) with eltype ComplexF64:
1.0 + 1.0im
2.0 - 1.0im

julia> reinterpret(reshape, Float64, s) # The actual memory is `([1.0, 2.0], [3.0, -1.0])`
2×2 reinterpret(reshape, Float64, StructArray(::Vector{Float64}, ::Vector{Float64})) with eltype Float64:
1.0 2.0
3.0 -1.0
```

If you already have a `StructArray`, the easiest way is to get the higher-dimensional
"struct-of-arrays" layout is to directly stack the components in memory order:

```julia
julia> using StackViews # lazily cat/stack arrays in a new tailing dimension

julia> StackView(StructArrays.components(s)...)
2×2 StackView{Float64, 2, 2, Tuple{Vector{Float64}, Vector{Float64}}}:
1.0 3.0
2.0 -1.0
```

StructArrays also provides a way to reconstruct from a given memory block via `dims` keyword:

```julia
julia> v = Float64[1 3; 2 -1]
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I use v = [1 3; 2 -1] then I get

julia> StructArray{ComplexF64}(v, dims=1)
0-element StructArray(StructArray(), StructArray()) with eltype ComplexF64 with indices 1:0

Is this expected?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you sure? I get

julia> v = Float64[1 3; 2 -1]
2×2 Matrix{Float64}:
 1.0   3.0
 2.0  -1.0

julia> StructArray{ComplexF64}(v, dims=1)
2-element StructArray(view(::Matrix{Float64}, 1, :), view(::Matrix{Float64}, 2, :)) with eltype ComplexF64:
 1.0 + 2.0im
 3.0 - 1.0im

which is the expected behavior. Btw, dims=ndims(v) would be the way to get contiguous views as component arrays (selecting on the last dimension).

Copy link
Member Author

@johnnychen94 johnnychen94 Aug 26, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry , I mean:

julia> v = Int[1 3; 2 -1]
2×2 Matrix{Int64}:
 1   3
 2  -1

julia> StructArray{ComplexF64}(v, dims=1)
0-element StructArray(StructArray(), StructArray()) with eltype ComplexF64 with indices 1:0

julia> StructArray{ComplexF64}(reinterpret(Float64, v), dims=1)
2-element StructArray(view(reinterpret(Float64, ::Matrix{Int64}), 1, :), view(reinterpret(Float64, ::Matrix{Int64}), 2, :)) with eltype ComplexF64:
 5.0e-324 + 1.0e-323im
 1.5e-323 + NaN*im

julia> v = ComplexF64[1 3; 2 -1]
2×2 Matrix{ComplexF64}:
 1.0+0.0im   3.0+0.0im
 2.0+0.0im  -1.0+0.0im

julia> StructArray{ComplexF64}(v, dims=1)
2-element view(::Matrix{ComplexF64}, 1, :) with eltype ComplexF64:
 1.0 + 0.0im
 3.0 + 0.0im

julia> StructArray{ComplexF64}(v, dims=2)
2-element view(::Matrix{ComplexF64}, :, 1) with eltype ComplexF64:
 1.0 + 0.0im
 2.0 + 0.0im

Interpreting the output is quite, hmmm, unintuitive. Maybe I just hit some undefined behaviors.

Copy link
Collaborator

@piever piever Aug 26, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, I see, it behaves a bit funny if the types don't match. It's probably because it has some logic to support nested cases, I've opened #200 to track this. Note that this constructor does not allocate, it can't use a matrix of integers to store components of ComplexF64.

2×2 Matrix{Float64}:
1.0 3.0
2.0 -1.0

julia> StructArray{ComplexF64}(v, dims=1) # the actual memory is `([1.0, 3.0], [2.0, -1.0])`
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, here and below the memory is always v, this constructor does not allocate.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh thanks for pointing this out! I didn't realize this, I only did a quick check @btime StructArray{ComplexF64}(v, dims=1) with small v = Float64[1 3; 2 -1] and I thought the memory allocations means copy 😂

julia> @btime StructArray{ComplexF64}(v, dims=1);
521.754 ns (8 allocations: 352 bytes)

2-element StructArray(view(::Matrix{Float64}, 1, :), view(::Matrix{Float64}, 2, :)) with eltype ComplexF64:
1.0 + 2.0im
3.0 - 1.0im

julia> s = StructArray{ComplexF64}(v, dims=2) # the actual memory is `([1.0, 2.0], [3.0, -1.0])`
2-element StructArray(view(::Matrix{Float64}, :, 1), view(::Matrix{Float64}, :, 2)) with eltype ComplexF64:
1.0 + 3.0im
2.0 - 1.0im
```

This, however, depends on the underlying data layout and how you interpret the memory block. You
should use this with caution because otherwise it might give you unexpected results.