Reading data from HDF5 file to StaticArray

I’m trying to read large datasets from HDF5 files into StaticArrays but I’m really struggling to understand how to get the data into a StaticArray.

I’m using the following function to read the HDF5 file and it works fine, returning the expected data:

function get_property(filepath::String, property::String)
    timestep_key = filename_base(filepath)
    fid = h5open(filepath, "r")
    timestep_key = keys(fid["TimestepData"])[1]
    query = read(fid["TimestepData"][timestep_key][property])
    # for testing purposes query is a (3, 44843) Matrix of Float64

    # Some code here that converts data to StaticArray?

    return query

Typically I always know one dimension (I’m getting force data so there would always be 3 components, so (3,n) size) so I thought it should be straightforward to convert to a StaticArray to take advantage of the faster computation speeds.

But I’m not able to do that conversion of the data before returning from the function.

I’ve tried:

> Internal Error: MethodError: no method matching *(::Type{StaticArrays.SVector{3, Float64}}, ::Matrix{Float64})

> Internal Error: DimensionMismatch("expected input array of length 3, got length 134529")

>Internal Error: MethodError: Cannot `convert` an object of type Matrix{Float64} to an object of type Float64
Closest candidates are:
  convert(::Type{T}, !Matched::Static.StaticFloat64{N}) where {N, T<:AbstractFloat} at C:\Users\jpmor\.julia\packages\Static\pkxBE\src\float.jl:26
  convert(::Type{T}, !Matched::Base.TwicePrecision) where T<:Number at C:\Users\jpmor\.julia\juliaup\julia-1.7.1+0~x64\share\julia\base\twiceprecision.jl:262
  convert(::Type{T}, !Matched::AbstractChar) where T<:Number at C:\Users\jpmor\.julia\juliaup\julia-1.7.1+0~x64\share\julia\base\char.jl:185

SMatrix{3,44843}(query) # This is still running after almost 1hr....

Can someone point out to me how to achieve this - I am assuming it is possible and simple, but the StaticArrays documentation is very basic and doesn’t seem to cover this.

If it’s a slow and computationally intensive task, then I will probably stick to a standard matrix, but I can’t get this to work at all.

A matrix of size (3, 44843) is much too large for StaticArrays. The StaticArrays README says:

Note that in the current implementation, working with large StaticArray s puts a lot of stress on the compiler, and becomes slower than Base.Array as the size increases. A very rough rule of thumb is that you should consider using a normal Array for arrays larger than 100 elements.

I’m not sure what query is supposed to represent, but if it is a collection of x,y,z points, or something else along the columns that you want to use together, you could consider using a Vector of SVector:

julia> using StaticArrays

julia> query = rand(3, 44843)
3×44843 Matrix{Float64}:
 0.62791    0.376845  0.0319929  …  0.254759  0.748335
 0.392331   0.469574  0.622981      0.154285  0.673124
 0.0994297  0.132254  0.791781      0.251181  0.802184

julia> [SVector{3,Float64}(col) for col in eachcol(query)]
44843-element Vector{SVector{3, Float64}}:
 [0.6279103058302006, 0.3923314542013667, 0.09942971962927671]
 [0.3768454590435031, 0.469573741906429, 0.13225405391447298]
1 Like

Thanks. I completely missed that statement about array size. It’s basically a matrix of forces (x,y,z) at nodes. There would be another set of points of (x,y,z) to provide the location of each node.

I had been thinking I could implement something like the FieldVector example in the StaticArray documentation for forces and positions.

Anyway, thanks!

If you want to use FieldVector just so you can access the x value with p.x, note that since StaticArrays 1.3, released 4 days ago, x/y/z/w field access is implemented for SVector and MVector as well:

julia> using StaticArrays

julia> p = SA[1.1, 2.2, 3.3]
3-element SVector{3, Float64} with indices SOneTo(3):

julia> p.x

julia> p.y

julia> p.z

For more info, see Named access with .{x,y,z,w} to elements of SVector and MVector by c42f · Pull Request #980 · JuliaArrays/StaticArrays.jl · GitHub.

1 Like

Oh that is nice. Something for free!

If you read in a huge array and don’t want to create a copy with SVectors, you can also reinterpret the original array like reinterpret(S, arr) where S is your chosen point representation.

1 Like

Sorry for asking, wouldn’t this other post be a solution for that?

Yes that looks like the same as what Julius proposed, only that example is with different dimensions, an SMatrix from a 3D Array. Here it is written out for this case:

julia> using StaticArrays

julia> query = rand(3, 44843)
3×44843 Matrix{Float64}:
 0.779144  0.664259  …  0.51808   0.953658
 0.535466  0.933785     0.767006  0.0611424
 0.337313  0.671113     0.293518  0.892914

julia> reinterpret(SVector{3, Float64}, query)
1×44843 reinterpret(SVector{3, Float64}, ::Matrix{Float64}):
 [0.779144, 0.535466, 0.337313]  …  [0.953658, 0.0611424, 0.892914]
1 Like