How to use differently sized static arrays per project

ufechner7 · June 13, 2021, 1:26pm

Hello,

I have a module called Utils (see KiteViewer/Utils.jl at sim · ufechner7/KiteViewer · GitHub).

It provides a function se() that returns a struct of simulation settings. These
settings are loaded from a YAML file.

Now I want to use different settings depending on the system that I am simulating.
But the project to use must be known BEFORE I load the module Utils.

The reason is that the settings contain the constant “segments” which defines
the size of static vectors I use for the simulation.

My current plan is to use a file project.txt witch contains the name of the
settings file that shall be loaded.

And if I change the settings during runtime the program shall update this file and
then it would re-start itself.

Is there a better way to achieve the desired behavior:

using statically sized arrays for the simulation
reload the program when the project and thus the size changes?

I a bit afraid that restarting might be pretty slow because it might cause a lot of recompilation.

Uwe

lmiq · June 13, 2021, 2:40pm

I think it is better to parameterize the structs with the size of the vectors. For example, this one:

struct ExtSysState
    time::Float64                          # time since launch in seconds
    orient::UnitQuaternion{Float32}        # orientation of the kite
    X::MVector{se().segments+1, MyFloat} 
...
end

Use instead:

struct ExtSysState{N,T}
    time::Float64                          # time since launch in seconds
    orient::UnitQuaternion{Float32}        # orientation of the kite
    X::MVector{N,T} 
...
end

then the definition of the vector size and type will occur in the initialization of the structure, not on module loading. If you follow this pattern in your code, the type of problem will be defined at run time each function.

ufechner7 · June 13, 2021, 6:09pm

I more-or-less implemented your suggestion, and it works.

See: KiteViewer/Utils.jl at sim · ufechner7/KiteViewer · GitHub

On the other hand the startup time of my app increased from 11s to 12s.

What was surprising for me was that I cannot do any calculations with the type parameter. I first tried to pass the number of segments, S and use S+1 as vector size, but that did not work. Now I use P, the number of points as type parameter.

lmiq · June 13, 2021, 6:32pm

What you mean here exactly? Three approaches could be used, I think:

Annotate the parametric types on the functions (this one is the most elegant, in my opinion, since the function will specialize to that type of a anyway:

function f(a::A{N,T}) where {N,T}
   # use N here
end

Get the size internaly:

function f(a)
   N = length(a.x) # a.x is the MVector{N,T}
   ...
end

Annotate he size in your struct:

struct A{N,T}
   N::Int
   x::Vector{N,T}
end

and initialize it accordingly (and use a.N there)

lmiq · June 13, 2021, 6:42pm

Some other thing: how many particles are you talking about? MVectors with large N (P in your case) are not a good idea and will make compile times be very large. Better data structures are standard vectors or, probably, a vector of SVectors.

ufechner7 · June 13, 2021, 6:44pm

P is currently 7, might be up to 14, but probably never larger…
I am pretty happy with these data structures for now, my inner loop needs about 500ns, with normal vectors it is 15 times as much.
You can find my inner loop in line 270ff of this file: KiteViewer/KPS3.jl at sim · ufechner7/KiteViewer · GitHub

For now it uses a constant for the number of segments, because if I do not use a constant it allocates… But I will experiment with type parameters, too, lets see…

ufechner7 · June 13, 2021, 6:55pm

This is the difficult code:

function residual!(res, yd, y, p, time)
    # unpack the vectors y and yd
    part = reshape(SVector{6*(SEGMENTS)}(y),  Size(3, SEGMENTS, 2))
    partd = reshape(SVector{6*(SEGMENTS)}(yd),  Size(3, SEGMENTS, 2))
    pos1, vel1 = part[:,:,1], part[:,:,2]
    pos = SVector{SEGMENTS+1}(if i==1 SVector(0.0,0,0) else SVector(pos1[:,i-1]) end for i in 1:SEGMENTS+1)
    vel = SVector{SEGMENTS+1}(if i==1 SVector(0.0,0,0) else SVector(vel1[:,i-1]) end for i in 1:SEGMENTS+1)
    posd1, veld1 = partd[:,:,1], partd[:,:,2]
    posd = SVector{SEGMENTS+1}(if i==1 SVector(0.0,0,0) else SVector(posd1[:,i-1]) end for i in 1:SEGMENTS+1)
    veld = SVector{SEGMENTS+1}(if i==1 SVector(0.0,0,0) else SVector(veld1[:,i-1]) end for i in 1:SEGMENTS+1)

If I do not define SEGMENTS as constant, it allocates.

lmiq · June 13, 2021, 7:28pm

Exactly. If one the input parameters contain the size as a type parameter, use something like:

and then it may become non allocating. The point is to make the function specialize to a given N, and then these sizes are known at compile time, even if not apparently constant for the user.

(“particles” in my field are usually on the thousands, so that is why I asked. But with those sizes that is fine)

Ps. For example: Question about allocations and generators - #3 by lmiq

ufechner7 · June 13, 2021, 7:39pm

Well, the function is a callback function of the IDA integrator, and I cannot change the function signature… But p is a parameter that I can freely choose (well, within curtain limits). So perhaps your suggestion will work…

ufechner7 · June 13, 2021, 8:07pm

So I tried:

function residual!(res, yd, y, p{N}, time) where N

and get the error message:

syntax: "p{N}" is not a valid function argument

So what exactly did you suggest?

Found it myself. This might work (at least no syntax error):

function residual!(res, yd, y, p::Vector{N}, time) where N

ufechner7 · June 13, 2021, 9:02pm

Final solution:

function residual!(res, yd, y::MVector{S, Float64}, p, time) where S
    # unpack the vectors y and yd
    part = reshape(SVector{S}(y),  Size(3, div(S,6), 2))
    partd = reshape(SVector{S}(yd),  Size(3, div(S,6), 2))
    pos1, vel1 = part[:,:,1], part[:,:,2]
    pos = SVector{div(S,6)+1}(if i==1 SVector(0.0,0,0) else SVector(pos1[:,i-1]) end for i in 1:div(S,6)+1)
    vel = SVector{div(S,6)+1}(if i==1 SVector(0.0,0,0) else SVector(vel1[:,i-1]) end for i in 1:div(S,6)+1)
    posd1, veld1 = partd[:,:,1], partd[:,:,2]
    posd = SVector{div(S,6)+1}(if i==1 SVector(0.0,0,0) else SVector(posd1[:,i-1]) end for i in 1:div(S,6)+1)
    veld = SVector{div(S,6)+1}(if i==1 SVector(0.0,0,0) else SVector(veld1[:,i-1]) end for i in 1:div(S,6)+1)

I could not use the p parameter, because I can use a float or a Vector of float for it, but not an SVector.

lmiq · June 13, 2021, 10:56pm

Sorry for that. It was a mistake because of lack os attention. What you finally did is correct, and I fixed the example to avoid confusion.