How to access columns from Vector{SVector{4,Float64}}?

General Usage
,
1

I have this script file. I originally planned to utilize my own ODE solver called RKF45 but convergence issues prevented that and I fell back on using ODEProblem and solve from DifferentialEquations. I had written my ODE function to utilize SVector, as that is what RKF45 uses, and I needed to do some rewriting to get my ODE function to work with DifferentialEquations. From solve, I get t and u where t is the t value vector and u is a solution Vector{SVector{4,Float64}} with 924 entries (corresponding to different values of t, the 4 elements of the SVector correspond to the 4 variables being integrated). My problem with this output is that I cannot seem to access individual columns, only individual elements or rows. The columns correspond to different integrated variables and I want to be able to access them so that I can plot them.

I have tried using SMatrix to convert u to a 924x4 SMatrix with Float64 elements with the code:

M = SMatrix{924,4,Float64}(u);

But this gave the error ERROR: DimensionMismatch: expected input array of length 3696, got length 924.

2

I have found that this code extracts the columns:

z = zeros(length(u), 1);
dz = zeros(length(u), 1);
theta = zeros(length(u), 1);
dtheta = zeros(length(u), 1);
for i in 1:length(u)
    ui = u[i];
    z[i]      = ui[1];
    dz[i]     = ui[2];
    theta[i]  = ui[3];
    dtheta[i] = ui[4];
end

But I am wondering if there’s a more efficient way to do so.

3

You could use

julia> using StaticArrays

julia> v = rand(SVector{4, Float64}, 924); v[1:3]
3-element Vector{SVector{4, Float64}}:
 [0.26989564837338975, 0.7293842397173794, 0.3764046148962866, 0.19753458799402568]
 [0.3685980243939463, 0.981358666933938, 0.242629242942751, 0.869518836692788]
 [0.33657278759867315, 0.4027440541365107, 0.45841596629600323, 0.15337648663757097]

julia> x = transpose(reshape(reinterpret(Float64, v), 4, :)); x[1:3, :]
3×4 Matrix{Float64}:
 0.269896  0.729384  0.376405  0.197535
 0.368598  0.981359  0.242629  0.869519
 0.336573  0.402744  0.458416  0.153376

julia> z = x[:, 1]; z[1:3]
3-element Vector{Float64}:
 0.26989564837338975
 0.3685980243939463
 0.33657278759867315

Here we

  • reinterpret to ‘forget’ we were grouping the Float64s in tuples of 4
  • reshape the resulting Vector{Float64} of length 3696 = 924 \cdot 4 into a Matrix. Because the four floats in an SVector{4, Float64} are next to each other in memory, and Julia is column major, we need a 4 \times 924 matrix
  • transpose to get an 924 \times 4 matrix.
4

You’re calling getindex per element, you can also do that with broadcasting:

julia> u = fill(SVector{4,Float64}(1:4), 3)
3-element Vector{SVector{4, Float64}}:
 [1.0, 2.0, 3.0, 4.0]
 [1.0, 2.0, 3.0, 4.0]
 [1.0, 2.0, 3.0, 4.0]

julia> z = getindex.(u, 1)
3-element Vector{Float64}:
 1.0
 1.0
 1.0

julia> dz = getindex.(u, 2)
3-element Vector{Float64}:
 2.0
 2.0
 2.0
2 Likes
5 
transpose(reshape(reinterpret(Float64, u)))

returns ERROR: DimensionMismatch: parent has 113648 elements, which is incompatible with size ().

6

HybridArrays.jl is nice for this.

7

You’re not supplying a dims(...) argument in reshape for the new size, which then seems to default to an empty tuple. Here you would need (4, length(u)) (or (4, :)).

1 Like
8

Note that this can be written more concisely as

reinterpret(reshape, Float64, v)'
3 Likes