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

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

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]

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))))

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

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?

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]

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)
3 Likes