Reinterpret examples for converting Julia multi-dimensional arrays and StaticArrays

I’m having trouble discerning how to convert between different representations of array structures using reinterpret. I’ve figured some things out, which I have posted below in case anyone else will find the examples useful. I end with a case that I am currently unable to solve and would appreciate any suggestions.

Suppose I have a type representing the Homogeneous coordinates of a planar point,

struct HomogeneousPoint{T <: AbstractFloat,N} 
    coords::NTuple{N, T}
end

and I instantiate it with 8 concrete points as follows:

pts2D_hom = map(HomogeneousPoint,
        [(4.244432486132958, 4.244432486132958, 1.0),
         (4.244432486132958, 16.755567513867042, 1.0),
         (8.244432486132958, 8.244432486132958, 1.0),
         (8.244432486132958, 12.755567513867042, 1.0),
         (12.755567513867042, 8.244432486132958, 1.0),
         (12.755567513867042, 12.755567513867042, 1.0),
         (16.755567513867042, 4.244432486132958, 1.0),
         (16.755567513867042, 16.755567513867042, 1.0)]).

This yields an Array of HomogeneousPoint

8-element Array{HomogeneousPoint{Float64,3},1}:
 HomogeneousPoint{Float64,3}((4.24443, 4.24443, 1.0))
 HomogeneousPoint{Float64,3}((4.24443, 16.7556, 1.0))
 HomogeneousPoint{Float64,3}((8.24443, 8.24443, 1.0))
 HomogeneousPoint{Float64,3}((8.24443, 12.7556, 1.0))
 HomogeneousPoint{Float64,3}((12.7556, 8.24443, 1.0))
 HomogeneousPoint{Float64,3}((12.7556, 12.7556, 1.0))
 HomogeneousPoint{Float64,3}((16.7556, 4.24443, 1.0))
 HomogeneousPoint{Float64,3}((16.7556, 16.7556, 1.0))

We can reinterpret this array as an 8 x 3 two-dimensional array with the command

mat = reinterpret(Float64,pts2D_hom ,(8,3))

resulting in

8×3 Array{Float64,2}:
  4.24443   1.0      12.7556 
  4.24443   8.24443   1.0    
  1.0      12.7556   16.7556 
  4.24443   1.0       4.24443
 16.7556   12.7556    1.0    
  1.0       8.24443  16.7556 
  8.24443   1.0      16.7556 
  8.24443  12.7556    1.0    

The rational for the construction, as oatlzzvztd points out in Reinterpret SVector, is that

the first argument to reinterpret should be the type of the element of the resulting array (i.e., the result of eltype(y).

Now suppose that I wanted to reinterpret pts2D_hom as an SMatrix instead. The closest I came up with is

smat = reinterpret(SMatrix{8,3,Float64,24},pts2D_hom,(1,))

which creates an array with a single entry of an SMatrix as follows

1-element Array{StaticArrays.SArray{Tuple{8,3},Float64,2,24},1}:
 [4.24443 1.0 12.7556; 4.24443 8.24443 1.0; … ; 8.24443 1.0 16.7556; 8.24443 12.7556 1.0]

However, what I am I to do if I don’t want an array containing an SMatrix, but I want the SMatrix on its own?

As far as I know, extracting it via getindex is the only way.
Note that arrays are heap allocated, while SMatrices are stack allocated.

You won’t change where the memory is by reinterpreting; you’re just changing the interpretation of the memory.
To get a fast stack-allocated matrix, you’ll have to actually move that data.

Of course, it’s probably just easier to do that the natural way

julia> x = randn(8,3);

julia> smx = SMatrix{8,3}(x)
8×3 SArray{Tuple{8,3},Float64,2,24}:
 -0.285864  -1.49749      0.232151
 -2.02672    0.213289     1.8493  
  0.550775  -0.565527     1.32064 
 -0.35845    0.00888187   1.15282 
  1.31989    1.17796      0.512057
  0.956868  -2.17719     -1.26094 
  1.64638    0.869573    -0.202741
 -0.209051  -1.15422     -0.991324

julia> @benchmark SMatrix{8,3}($x) #Does not allocate (heap) memory
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     37.959 ns (0.00% GC)
  median time:      37.973 ns (0.00% GC)
  mean time:        39.665 ns (0.00% GC)
  maximum time:     89.337 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     992

julia> function smreinterpret(x, ::Val{nrow}, ::Val{ncol}) where {ncol, nrow}
           @inbounds out = reinterpret(SMatrix{nrow,ncol,Float64,nrow*ncol},vec(x),)[1]
           out
       end
smreinterpret (generic function with 1 method)

julia> smreinterpret(x, Val{8}(), Val{3}())
8×3 SArray{Tuple{8,3},Float64,2,24}:
 -0.285864  -1.49749      0.232151
 -2.02672    0.213289     1.8493  
  0.550775  -0.565527     1.32064 
 -0.35845    0.00888187   1.15282 
  1.31989    1.17796      0.512057
  0.956868  -2.17719     -1.26094 
  1.64638    0.869573    -0.202741
 -0.209051  -1.15422     -0.991324

julia> @benchmark smreinterpret($x, Val{8}(), Val{3}())
BenchmarkTools.Trial: 
  memory estimate:  80 bytes
  allocs estimate:  2
  --------------
  minimum time:     409.875 ns (0.00% GC)
  median time:      436.460 ns (0.00% GC)
  mean time:        479.329 ns (8.49% GC)
  maximum time:     285.814 μs (99.83% GC)
  --------------
  samples:          10000
  evals/sample:     200

1 Like