Reinterpret array into SMatrices of different sizes for Interpolations. Julia v0.6 vs now

General Usage
1

Hi all,
I’m trying to check whether my approach to translation something that used to work in previous julia versions is the best way to go.

I have the following use case. I have values stored in a 4-dimensional array y. I want to leverage this awesome functionality of Interpolations.jl whereby one can interpolate multiple functions over the same domain(s) in a single call. I want to interpolate dimensions 3 and 4 of x, obtaining a 2D SMatrix as a result each time (basically the interpolated slices of dimensions 1 and 2 of y). Here is the example that works in Julia 0.6. (real application uses Float, but that’s not important):

# julia 0.6
julia> y = reshape(collect(1:48),2,2,3,4)

# that command here does not work in v1.9 - see below
julia> z = reinterpret(SMatrix{2,2,Int,2*2}, y, (3,4))
3×4 Array{StaticArrays.SArray{Tuple{2,2},Int64,2,4},2}:
 [1 3; 2 4]     [13 15; 14 16]  [25 27; 26 28]  [37 39; 38 40]
 [5 7; 6 8]     [17 19; 18 20]  [29 31; 30 32]  [41 43; 42 44]
 [9 11; 10 12]  [21 23; 22 24]  [33 35; 34 36]  [45 47; 46 48]

julia> i = interpolate(([1,2,3],[1,2,3,4]), z, Gridded(Constant()))
3×4 interpolate((::Array{Int64,1},::Array{Int64,1}), ::Array{StaticArrays.SArray{Tuple{2,2},Int64,2,4},2}, Gridded(Constant())) with element type StaticArrays.SArray{Tuple{2,2},Float64,2,4}:
 [1 3; 2 4]     [13 15; 14 16]  [25 27; 26 28]  [37 39; 38 40]
 [5 7; 6 8]     [17 19; 18 20]  [29 31; 30 32]  [41 43; 42 44]
 [9 11; 10 12]  [21 23; 22 24]  [33 35; 34 36]  [45 47; 46 48]

julia> i[1,1]
2×2 StaticArrays.SArray{Tuple{2,2},Int64,2,4}:
 1  3
 2  4

julia> i[2,3]
2×2 StaticArrays.SArray{Tuple{2,2},Int64,2,4}:
 29  31
 30  32

so far, so good. now, reinterpret has changed it seems. Here is my approach in modern Julia, please advise whether that’s ok. thanks.

# julia 1.9
julia> y = reshape(collect(1:48),2,2,3,4);
julia> z = reinterpret(SMatrix{2,2,Int,4},reshape(y,12,4))  # notice the reshape!
3×4 reinterpret(SMatrix{2, 2, Int64, 4}, ::Matrix{Int64}):
 [1 3; 2 4]     [13 15; 14 16]  [25 27; 26 28]  [37 39; 38 40]
 [5 7; 6 8]     [17 19; 18 20]  [29 31; 30 32]  [41 43; 42 44]
 [9 11; 10 12]  [21 23; 22 24]  [33 35; 34 36]  [45 47; 46 48]

julia> i = interpolate(([1,2,3],[1,2,3,4]), z, Gridded(Constant()))
3×4 interpolate((::Vector{Int64},::Vector{Int64}), ::Matrix{SMatrix{2, 2, Int64, 4}}, Gridded(Constant{Nearest}())) with element type SMatrix{2, 2, Float64, 4}:
 [1 3; 2 4]     [13 15; 14 16]  [25 27; 26 28]  [37 39; 38 40]
 [5 7; 6 8]     [17 19; 18 20]  [29 31; 30 32]  [41 43; 42 44]
 [9 11; 10 12]  [21 23; 22 24]  [33 35; 34 36]  [45 47; 46 48]

julia> i(1,1)
2×2 SMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
 1  3
 2  4

julia> i(2,3)
2×2 SMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
 29  31
 30  32
1 Like
2

I think this is correct. Another way of writing this is

julia> z = reinterpret(reshape, SMatrix{2,2,Int,4}, reshape(y, :, 3, 4))
3×4 reinterpret(reshape, SMatrix{2, 2, Int64, 4}, ::Array{Int64, 3}) with eltype SMatrix{2, 2, Int64, 4}:
 [1 3; 2 4]     [13 15; 14 16]  [25 27; 26 28]  [37 39; 38 40]
 [5 7; 6 8]     [17 19; 18 20]  [29 31; 30 32]  [41 43; 42 44]
 [9 11; 10 12]  [21 23; 22 24]  [33 35; 34 36]  [45 47; 46 48]
3

If the use of the reinterpret function is not mandatory, there are alternatives

x1=SMatrix{2,2,Int}.(eachslice(y;dims=(3,4)))
x2=[SMatrix{2,2,Int}(m) for m in (eachslice(y;dims=(3,4)))]
x3=map(m->SMatrix{2,2,Int}(m), (eachslice(y;dims=(3,4))))
4

Oh yes -I never understood the use if the first reshape there. Thanks!

5

Interesting! No, reinterpret is not necessary, just haven’t thought about the others. Nice. But I guess that is allocating whereas reinterpret maps the same memory in the original array?

6

more or less, everyone allocates

julia> @btime x3=map(m->SMatrix{2,2,Int}(m), (eachslice(y;dims=(3,4))))
  450.761 ns (7 allocations: 704 bytes)
3×4 Matrix{SMatrix{2, 2, Int64, 4}}:
 [1 3; 2 4]     [13 15; 14 16]  [25 27; 26 28]  [37 39; 38 40]
 [5 7; 6 8]     [17 19; 18 20]  [29 31; 30 32]  [41 43; 42 44]
 [9 11; 10 12]  [21 23; 22 24]  [33 35; 34 36]  [45 47; 46 48]

julia> @btime convert.(SMatrix{2,2,Int},eachslice(y;dims=(3,4)))
  398.000 ns (7 allocations: 704 bytes)
3×4 Matrix{SMatrix{2, 2, Int64, 4}}:
 [1 3; 2 4]     [13 15; 14 16]  [25 27; 26 28]  [37 39; 38 40]
 [5 7; 6 8]     [17 19; 18 20]  [29 31; 30 32]  [41 43; 42 44]
 [9 11; 10 12]  [21 23; 22 24]  [33 35; 34 36]  [45 47; 46 48]

julia> @btime reinterpret(reshape, SMatrix{2,2,Int,4}, reshape(y, :, 3, 4))
  240.769 ns (4 allocations: 160 bytes)
3×4 reinterpret(reshape, SMatrix{2, 2, Int64, 4}, ::Array{Int64, 3}) with eltype SMatrix{2, 2, Int64, 4}:
 [1 3; 2 4]     [13 15; 14 16]  [25 27; 26 28]  [37 39; 38 40]
 [5 7; 6 8]     [17 19; 18 20]  [29 31; 30 32]  [41 43; 42 44]
 [9 11; 10 12]  [21 23; 22 24]  [33 35; 34 36]  [45 47; 46 48]


julia> @btime reinterpret(reshape, SMatrix{2,2,Int,4}, reshape($y, :, 3, 4))
  27.035 ns (2 allocations: 96 bytes)
3×4 reinterpret(reshape, SMatrix{2, 2, Int64, 4}, ::Array{Int64, 3}) with eltype SMatrix{2, 2, Int64, 4}:
 [1 3; 2 4]     [13 15; 14 16]  [25 27; 26 28]  [37 39; 38 40]
 [5 7; 6 8]     [17 19; 18 20]  [29 31; 30 32]  [41 43; 42 44]
 [9 11; 10 12]  [21 23; 22 24]  [33 35; 34 36]  [45 47; 46 48]

julia> @btime convert.(SMatrix{2,2,Int},eachslice($y;dims=(3,4)))
  71.942 ns (1 allocation: 448 bytes)
3×4 Matrix{SMatrix{2, 2, Int64, 4}}:
 [1 3; 2 4]     [13 15; 14 16]  [25 27; 26 28]  [37 39; 38 40]
 [5 7; 6 8]     [17 19; 18 20]  [29 31; 30 32]  [41 43; 42 44]
 [9 11; 10 12]  [21 23; 22 24]  [33 35; 34 36]  [45 47; 46 48]
7

The allocations in reinterpret are minimal, arising from the array reshape. The allocations in the broadcast/map operations are of the order of the array. In this example, the small size of the array is masking the impact. If we choose larger trailing sizes:

julia> y = reshape(collect(1:prod((2,2,30,40))),2,2,30,40);

julia> @btime convert.(SMatrix{2,2,Int},eachslice($y;dims=(3,4)));
  3.408 μs (2 allocations: 37.55 KiB)

julia> @btime map(m->SMatrix{2,2,Int}(m), (eachslice($y;dims=(3,4))));
  5.986 μs (2 allocations: 37.55 KiB)

julia> @btime reinterpret(reshape, SMatrix{2,2,Int,4}, reshape($y, :, size($y)[3:4]...));
  35.289 ns (2 allocations: 96 bytes)
4 Likes