Agents.jl: Use Static Vector for position type of GridSpace and ContinuousSpace?

[Datseris]

We are going to 3.0, so maybe it is time to consider this seriously. Using tuples has drawbacks:

1. Basic linalg functions don't work on them, e.g. norm.
2. We can't use Distances.jl
3. Tuples do not extend the vector API, making users have difficulties.

The drawback in making the change is:

1. It will take a bunch of development time :D
2. Users will have to initialize static vectors for agent positions. Thankfully, the shorthand syntax SA[1, 2] makes a static vector. It is not that far from making a tuple!

[Datseris]

ah, a benefit of static arrays: x+y works without need to broadcast if x and y are the same type.

[Libbum]

The benefits here are mostly ease of future development, and performance (albeit in limited places), right? From the user's point of view SA[1, 2] is not that big of a deal vs (1,2). However, how often do we expect users to want to get mean(a.pos)? We use a.pos[1] or first(a.pos) in some examples for data collection, but I see these as contrived due to our test agents being so simple.

Other than the euclidean/cityblock metrics we've discussed earlier would make sense to include here? Is Distances.jl really needed? Do we want to provide the entire suite of exotic distance metrics?

What other areas do we see in the pro column for this change? Is it worth the overhead in the long term?

If we can answer those questions satisfactorily, this seems like a reasonable step forward, and I doubt it'll take anywhere near the amount of effort the data collection did.

[Datseris]

However, how often do we expect users to want to get mean(a.pos)?

I don't think it makes much sense to data collect mean(a.pos), but I think it is super useful to be able to collect norm(a.vel), or norm(a.pos - b.pos). In fact, in @itsdfish 's flocking example, it might even be useful to use the dot product dot(a.vel, b.vel), to calculate the angle betweeen the direction of two birds. Having access to vector operations can only help in the forseeable future.

Is Distances.jl really needed?

Not really, but it is more sustainable. If I could use distances in the source code, the code would be smaller. We are not the ones "providing" the exotic distances, Distances.jl is. Julia is a language that allows such powerful communication between packages, so we should leverage it. There are dedicated professionals that have made all these distance functions and tested them, we could re-use their effort instead.

[Datseris]

(actually, if I didn't have to work on SignalDecomposition.jl, I'd probably do this change today)

[Libbum]

Solid enough motivation (was being devil's advocate for the most part).

I'm looking at the structure of #113 and building Daisyworld currently, but not sure how much time I have today.

[itsdfish]

Just to follow up, many of the concerns and problems I initially encountered were resolved, at least for the flocking model. My initial concern was that tuples would create unnecessary allocations for velocity calculations. That turned out not to be the case, as @Datseris informed me. After reworking the model, it was fairly easy to accommodate tuples. However, there are still some cases where some linear algebra operations might not be possible. I looked at a few common operations below.

Perhaps one approach is to make the position and velocity slots more general and use tuples as a default, so that convenience functions can still be used. The user can override the default tuple type and perform some of the initial setup that is normally automated with convenience functions. One way I have done this in the past is to define a union type and use parametric types to make the position field concrete:

Pos = Union{NTuple,SVector}
mutable struct MyType{P<:Pos}
    pos::P
end

Here are some examples where tuples work and do not work.

x is a tuple and y is a vector in the following examples.

Outer products do not work with tuples

julia> y*y'
2×2 Array{Float64,2}:
 0.733248  0.232088
 0.232088  0.0734605

julia> x*x'
ERROR: MethodError: no method matching adjoint(::Tuple{Float64,Float64})
Closest candidates are:
  adjoint(::Missing) at missing.jl:100
  adjoint(::LightGraphs.DefaultDistance) at /home/dfish/.julia/packages/LightGraphs/UPjU9/src/distance.jl:22
  adjoint(::Number) at number.jl:169
  ...
Stacktrace:
 [1] top-level scope at none:0

dot product works with tuples if the dot function is called

julia> y'*y
0.8067087955687413

julia> dot(y,y)
0.8067087955687413

julia> dot(x,x)
0.18

julia> x'*x
ERROR: MethodError: no method matching adjoint(::Tuple{Float64,Float64})
Closest candidates are:
  adjoint(::Missing) at missing.jl:100
  adjoint(::LightGraphs.DefaultDistance) at /home/dfish/.julia/packages/LightGraphs/UPjU9/src/distance.jl:22
  adjoint(::Number) at number.jl:169
  ...
Stacktrace:
 [1] top-level scope at none:0

Contrary to what I initially reported, norm does work:

julia> norm(y)
0.8981696919673594

julia> norm(x)
0.4242640687119285

matrix multiplication does not work:

julia> m*x
ERROR: MethodError: no method matching *(::Array{Float64,2}, ::Tuple{Float64,Float64})
  ...
Stacktrace:
 [1] top-level scope at none:0

julia> m*y
2-element Array{Float64,1}:
 0.9849323911336143
 0.3020894832868109

[Datseris]

oh wow, actually _most_ operations do work fine. perhaps there is no reason to spend effort to do such a change if the only benefit is using distances (and we actually have already implemented the two most commonly used distances anyway).

Perhaps one approach is to make the position and velocity slots more general and use tuples as a default, so that convenience functions can still be used.

This is too complicated to be useful in my eyes and there is not much benefit in allowing such freedom.

---

I think we can close this issue, unless someone is really into implementing this change (I definitely don't have the time for it)

[itsdfish]

Yeah. 2 of the 4 common operations will work. Usually making code more generic is the preferred approach, but sometimes that is not feasible. I'm fine with closing the issue.

I think the next best alternative for operations like outer products and matrix multiplication would be to copy the velocity tuple to a vector and copy back to a tuple after the operations are performed.

[Datseris]

you can directly convert to and from static arrays:

julia> x = (0.5, 0.5)
(0.5, 0.5)

julia> SVector(x)
2-element SArray{Tuple{2},Float64,1,2} with indices SOneTo(2):
 0.5
 0.5

julia> Tuple(ans)
(0.5, 0.5)

[itsdfish]

Even better. It does not seem to allocate new memory. This might be good to place in a performance tips section.

[Datseris]

It does not seem to allocate new memory

We already discussed this on Slack: making immutable stuff doesn't allocate.

This might be good to place in a performance tips section.

I don't agree. This is general Julia knowledge, and we shouldn't try to teach this in Agents.jl, as they are entirely orthogonal. There would be no end to such tips, and in fact they would probably be as large as our entire tutorial. For example, have a look at https://docs.julialang.org/en/v1/manual/performance-tips/

[itsdfish]

I don't agree. This is general Julia knowledge, and we shouldn't try to teach this in Agents.jl, as they are entirely orthogonal.

I don't know if this is general Julia knowledge. Its ultimately an empirical question. However, as far as I can tell, the memory allocation behavior of immutables is not discussed in the performance section. I suspect this would not be clear to many users.

There would be no end to such tips, and in fact they would probably be as large as our entire tutorial.

Yeah. I agree that it would be cumbersome and unproductive to document a comprehensive list of performance tips. Just to be clear, that is not what I was suggesting. However, I think there might be a limited set of performance tips that would be helpful for users. The reason I suggested including this in a performance tips section is because models with a spatial or other physical components may need to use linear algebra operations that are not compatible with tuples. The most performant solution appears to be StaticArrays, which is an external package that users may not be aware of and is not suggested by the error message. As a result, this question might be asked again or a user might adopt a suboptimal workaround.

I think the use of Unions/Abstract containers for heterogeneous agents is another useful performance tip or observation. Although the Julia performance tips mentions that abstract containers cause performance decrements, it provides little information about the magnitude of the problem or what factors moderate the performance decrement. The benchmark I performed last weekend revealed that the slowdown was only a factor of 2-4 in fairly extreme cases and depends on the number of interactions between agents. In many cases, such as the sheep-wolf model, the implications for performance will be smaller and perhaps negligible.

[Datseris]

Yeah. I agree that it would be cumbersome and unproductive to document a comprehensive list of performance tips. Just to be clear, that is not what I was suggesting. However, I think there might be a limited set of performance tips that would be helpful for users. The reason I suggested including this in a performance tips section is because models with a spatial or other physical components may need to use linear algebra operations that are not compatible with tuples. The most performant solution appears to be StaticArrays, which is an external package that users may not be aware of and is not suggested by the error message. As a result, this question might be asked again or a user might adopt a suboptimal workaround.

If you can squeeze a tip as a single sentence in the docstring of ContinuousSpace, feel free to open a PR for it!

The benchmark I performed last weekend revealed that the slowdown was only a factor of 2-4 in fairly extreme cases and depends on the number of interactions between agents.

There is a reason for that: almost no operation is contaminated by the type instability besides "getting the agent out of the dictionary". This happens because of the huge amount of function barriers in the source code of Agents.jl. So, in a sense it is because of good design that you have such a minor impact. In other scenarios it could be much larger. Also, later Julia versions have optimizations for small Unions, which also helps our scenario here. And if you have a union of parametrically typed agents, it would be even worse.

[itsdfish]

Agreed. I think you guys did a good job of minimizing the impact of heterogeneous types. If you think the benchmark results are useful beyond developer knowledge, I would be happy to submit a short write up.

I'll submit a PR with a short reference to StaticArrays.