VectorOfArrays: iterating over pairs; type qualifier when containing specific SVector

New to Julia
, , ,
1

I’m working on a project that involves solving something like a multi-particle problem (the application is not terribly important for my question), where each particle moves under a velocity field that depends on locations of all particles in the system (think planets affecting each other by gravity). I’m building the velocity field function that would eventually get passed onto ODEProblem from DifferentialEquations.jl.

To store each particle (in real plane, at least for now), I use the alias to

const Vec2D = StaticArrays.SVector{2}

I’m trying to use VectorOfArray from RecursiveArrayTools.jl to represent the ensemble of particles. (UPDATE: I originally used a plain vector of static vectors but ran into trouble with adaptive ODE integrators. Based on ODE of Vector of Matrices I decided to try out VectorOfArray.)

I’m new convert from MATLAB to Julia, I’m adjusting to the multiple dispatch and types and syntax, so I have a few questions. (Any other suggestions based on the limited code I’m posting here are more than welcome.)

For example, I initialize two ensembles, of N and M random particles by:

pN = VectorOfArrays( [Vec2D(c) for c in eachcol(rand(2,N))] )
pM = VectorOfArrays( [Vec2D(c) for c in eachcol(rand(2,M))] )

Q1 I want to perform a computation over all pairs of particles. Guided by Broadcasting over combinations of inputs

For example, to simply list all pairs (without computing on them), I hoped that the following should work:

[ (v,w) for v in pN, w in pM ] # or
[ args for args in Iterators.product(pN, pM)]

but I get

ERROR: MethodError: no method matching Array{Tuple{SVector{2, Float64}, SVector{2, Float64}}, 4}(::UndefInitializer, ::Tuple{Int64, Int64})
Closest candidates are:
  Array{T, N}(::UndefInitializer, ::Tuple{Vararg{Int64, N}}) where {T, N} at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/boot.jl:469
  Array{T, N}(::UndefInitializer, ::Tuple{Vararg{Integer, N}}) where {T, N} at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/baseext.jl:36
  Array{T, N}(::Missing, ::Any...) where {T, N} at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/baseext.jl:43

So I’ve done

[ (v,w) for v in pN[:], w in pM[:] ] # or
[ args for args in Iterators.product(pN[:], pM[:])]

But using [:] seems a bit… clunky? Is this the best way to go about doing this, or can I avoid the VectorOfArray → Vector conversion ([:])? Should I be even worried (from efficiency/performance standpoint) about this?

Q2. Say that I wanted to define a type alias Ensemble to ensure it’s a VectorOfArrays holding only Vec2D, how would I do it?

My guesses

const Ensemble = VectorOfArrays{Vec2D}

or

const Ensemble{T} = VectorOfArray{T,2,Vector{Vec2D{T}} } where {T <: Number }

(second guided by typeof(pN) for pN defined as above) did not work. They evaluate ok, but when I try a drop-in replacement

pN = Ensemble( [Vec2D(c) for c in eachcol(rand(2,N))] )

I get a

ERROR: MethodError: no method matching (VectorOfArray{T, 2, Array{SVector{2, T}, 1}} where T<:Number)(::Vector{SVector{2, Float64}})

and similar for the first type.

2

Is there any specific reason for the use of VectorOfArrays? Why not simply a vector of static vectors?

(That would make all your first attempts work)

3

Thanks for the suggestion - that’s what I originally had in my code, but I’ve been having a hard time getting a plain vector of static vectors to play nice with adaptive ODE integrators. Guided by discussion ODE of Vector of Matrices I decided to try out VectorOfArray.

If I were doing this in MATLAB, I would simply store them as a matrix (each column one particle), but I’m trying to “stretch” myself and try out other ways that (at least to me) appear more julian.

4

I think this is the way to go, and you can reinterpret the vector of static vectors as one back and forth without cost.

I’m no user of the VectorOfArrays alternative, but it seems to me that it is thought for cases where the elements of the array are mutable objects. For very small vectors, where SVectors are the best representation, that doesn’t seem appropriate.

5

A Vector of StaticVectors is equivalent in memory representation to a matrix, so in theory it should work.

But yes, I agree it’s a bit clunky sometimes to use VectorOfArray because its indexing is A[i,j] but then A[i] is the iteration of the vectors. This can break some generic assumptions. It was an idea in 2016 which I don’t think panned out: I think we need to change it so its indexing is always scalar. Though that would break the interface in DifferentialEquations.jl so it’ll need a major version bump there. We will need to do that sometime though, that and change of retcodes from Symbol to Enum are what are currently slated as the DiffEq v8 changes. Once that indexing change is done, all of your examples here should work without any workarounds.

In the meantime, your case might be better done with views of a Matrix, which also shouldn’t allocate. Or you can use reinterpret to change a Matrix into an array of staticvectors at any point without memory allocates (since it’s the same bit layout), so that’s another possibility.

1 Like
6

Thank you both! My previous solution involved @viewing the matrices, and I moved away from that, but I could go back. (I learned a few things about multiple dispatch and non-allocations along the way, so maybe I’ll do better than my first try.)

Would/could either of you answer my second question? Regardless of the VectorOfArray not being the best solution, could either of you (or someone else) say why the type aliases

const Ensemble = VectorOfArrays{Vec2D}

or

const Ensemble{T} = VectorOfArray{T,2,Vector{Vec2D{T}} } where {T <: Number }

don’t work as a drop in replacement in

pN = VectorOfArrays( [Vec2D(c) for c in eachcol(rand(2,N))] )

? I’m guessing I am missing something key about specifying parameters of types, but I couldn’t quite get it from the docs or permuting where I place the braces blindly.

7

I’m not sure, but maybe the problem is that SVector{2} is still abstract (in opposition to SVector{2,Float64} for example).