Agents.jl: Function to get all pairs within interaction radius

[lhupe]

@Datseris asked me to implement a simple model for growing dense bacterial colonies to test the usability of the new continuous space API.

In a dense bacterial colony (or indeed any dense granular system), a single agent can interact mechanically with several other agents at the same time. However, interacting_pairs only returns pairs of nearest neighbours. It would be nice to have a function that actually returns all pairs of agents within one interaction radius of each other.

[Datseris]

To be clear, space_neighbors is not good for your use case because for each interacting pair you only want to do the interaction once, correct? Having a look at the source code of interacting_pairs, it seems possible, to create a similar function that doesn't use nearest_neighbor but space_neighbors directly. In that source I have to keep track the pair IDs manually anyway, so it shouldn';t be that much harder.

[Datseris]

@lhupe how do you imagine the API of the function that you ask for here? I am having trouble to imagine what it should do and what it should return. Actually, I am not yet sure what are the input arguments.

[lhupe]

@Datseris I think there was some confusion on my part about the exact function of interacting_pairs!.
Does it return an iterator over all pairs of particles that are within one r neighborhood of each other or just over pairs of nearest neighbours?

I am getting some strange behaviour, which originally caused me to assume that the latter case is true. However, on reexamining my results, I find that the behaviour is inconsistent wih either hypothesis.

In the picture below, you can see a configuration of four of my Agents (plotted using a custom function). I have visualized the output of interacting_pairs! (with a sufficiently large interaction radius) on this configuration with the connecting lines. As you can see, the connection from the beige to the cyan agent is missing.

[Libbum]

Can you send us an MWE of this four agent cluster? We need to build some tests around this scenario I think. Then we can identify if interacting_pairs! needs to be modified (/fixed?) or we require a new method.

[Datseris]

Hi Lukas, thanks for showing this, it means that interacting_pairs is bugged. You should not have a connection between orange and magenta. Let's see the current source:

function interacting_pairs(model, r)
  pairs = Tuple{Int, Int}[]
  for id in keys(model.agents)
    # Skip already checked agents
    any(isequal(id), p[2] for p in pairs) && continue
    a1 = model[id]
    a2 = nearest_neighbor(a1, model, r)
    a2 ≠ nothing && push!(pairs, (id, a2.id))
  end
  return PairIterator(pairs, model.agents)
end

The PairIterator is not that relevant now, it just iterates over IDs but returns Agents instead of IDs because it is more convenient.

What I've tried to do with this function is to make it so each Agent belongs to only one pair, it;s nearest neighbor. Now, this is a bit ambiguous, because depending on if you started the process with the blue agent, you get a blue-orange and cyan-magenta pair. If you started with the orange agent, you'd get an orange-cyan and blue-magenta pair (always for sufficiently large r). The simple reason is that Orange's nearest neighbor is not Blue, but Blue's nearest neighbor is Orange.

I don't think there is a performant way to get rid of this inconsistency besides calculating the full NxN distance matrix (N = nagents(model)), and then sorting it by distance.

[Libbum]

There is an O(N) rather than O(N^2) method that I'm aware of called D-Cell. I haven't had the chance to finish my implementation of it in spherical-cow, but we can perhaps investigate the process at a later date. For now we could just consider the O(N^2) implementation and push a performance update in v3.x.

Key literature

[Datseris]

Yeah I can work on that maybe over the weekend if I find the time... And I agree that no reason to push for anything other than the absolute minimum. But I also want to bugfix the current implementation because it is the fastest one and sometimes one doesn't not necessarily care about this inaccuracy, depending on the dynamics of the system (e.g. it doesn't matter for where I used it, the social distancing example).

The thing is, even when we are done fixing this, it doesn't help what the user asks in the original question. And I don't know how to do this generically either.

@lhupe this is a way to solve your problem:

1. add a field interactedwith::Vector{Int} to your agents.
2. in agent_step! use space_neighbors. interact each agent with its space neighbors, and push in their fields interactedwith the original agent's id .
3. Loop over agents, finding space_neighbors but for each id in space neighbors check if it has already interacted. If yes, skip interaction. If no, do same as step 2.
4. After a the loop, add a model_step! function that simply clears all interactedwith fields.

This seems easy. I wonder if we could support such behavior directly from agents... @lhupe I would still appreciate you giving me some idea of how the API of the original question would look like. What would the function you want take as an input and what would it return?

In other news, it seems like the update_vel! functionality of ContinuousSpace doesn't help you here @lhupe ? When I coded that I had your system in mind but maybe it doesn't make sense, as the forces come from other agents. update_vel! would be helpful if e.g. the agents are magnetic and there is a perpendicular magnetic field.

[lhupe]

1. What I originally asked for is a function that – given the same input as interacting_pairs returns an iterator over all pairs of agents within one interaction radius of each other. In the example above, since the interaction radius is slightly larger than the yellow/magenta distance, it should return all the pairs highlighhted in the image as well as the yellow/cyan pair.
2. In fact, I could't even use the default move_agent!() function, because each agent is described by five dynamic variables (2D position, length, orientation, growth progress), so I wrote a custom agent_step!. However, this function is still seperated into calculating the "velocities" and actually evolving the parameters.

[Datseris]

Okay, I think (1) is ~easy~ possible to do based on the definition of interacting_pairs I cite above. For (2), I guess its because you want to update the velocities of all agents at the same time and then move all agents at the same time?

[Datseris]

@lhupe I can't reproduce your problem. Can you give specific agent coordinates? I do:

  model = ABM(Agent6, ContinuousSpace(2))
  for i in 1:4
    y = i == 2 ? 0.25 : 0.0
    add_agent_pos!(Agent6(i, ((i-1)*1.0, y), (0.0, 0.0), 0), model)
  end

julia> interacting_pairs(model, 1.0).pairs
2-element Array{Tuple{Int64,Int64},1}:
 (4, 3)
 (2, 1)

(and also clarify which metric you use)

[lhupe]

These are the agents positions, sorted by their IDs

(7.074386436066224, 4.963014649338054)
(5.831962448496828, 4.926297135685473)
(5.122087781793935, 5.300031210394806)
(3.9715633336430156, 4.8106570045816675)

This is the line I use to construct the world

space = ContinuousSpace(2, extend = (10, 10), periodic = false, metric = "euclidean")

The interaction radius for producing the above image was 2.0.

[Datseris]

I can't reproduce:

using Agents
mutable struct Agent6 <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    weight::Float64
end
  space = ContinuousSpace(2, extend = (10, 10), periodic = false, metric = :euclidean)
  model = ABM(Agent6, space)
  pos = [
    (7.074386436066224, 4.963014649338054)
    (5.831962448496828, 4.926297135685473)
    (5.122087781793935, 5.300031210394806)
    (3.9715633336430156, 4.8106570045816675)
  ]
  for i in 1:4
    add_agent_pos!(Agent6(i, pos[i], (0.0, 0.0), 0), model)
  end
  pairs = interacting_pairs(model, 2.0).pairs

julia> interacting_pairs(model, 2.0).pairs
2-element Array{Tuple{Int64,Int64},1}:
 (4, 3)
 (2, 1)

[Datseris]

Even if I run the pair collection algorithm with reverse IDs I still see no difference.

[lhupe]

I think I've found the key: the agents in the above example had IDs 3, 4, 5, 6 (because they were second-generation cells spawned by the divisision of cells 1 & 2 respectively). If you modify your MWE to

    add_agent_pos!(Agent6(i+2, pos[i], (0.0, 0.0), 0), model)

you should be able to reproduce the bug.

[Datseris]

Not really, doing the i+2 gives

julia> interacting_pairs(model, 2.0).pairs
2-element Array{Tuple{Int64,Int64},1}:
 (6, 5)
 (4, 3)

It also makes not much sense that increasing all IDs by 2 would change the behavior of the function. Looking at the source it is clear that only the operation == is used on IDs. So as long as there are no equal ids there is no problem.

Maybe you have created two agents with the same ID?

[lhupe]

Huh.
When I run the modified MWE, I get

3-element Array{Tuple{Int64,Int64},1}:
 (4, 3)
 (5, 3)
 (6, 5)

[Datseris]

God damn, you are right. I was testing many different schedulers, and turns out only the default one (fastest) does this. I'll look more into it and report.

[Libbum]

julia> pairs = interacting_pairs(model, 2.0).pairs
3-element Array{Tuple{Int64,Int64},1}:
 (4, 3)
 (5, 3)
 (6, 5)

I'm getting the same as @lhupe with the modification

[Datseris]

I've found the bug: The function does not check if an already found 2nd part of the pair happens to be a nearest neighbor of a new id.

I am pushing a PR for each agent belonging only to a single unique pair.

Thee function Lukas originally asked for is left as a TODO.

[Libbum]

@lhupe: can you check #225 and let me know if the behaviour you expect is the output of interacting_pairs(model, 2.0; method = :all)?

[Libbum]

Latest master has this functioning now.

interacting_pairs(model, 2.0) for nearest neighbor pairs within the radius,
interacting_pairs(model, 2.0; all = true) for all pairs within radius regardless of neighbor status.

[Libbum]

There were some issues with the previous implementation, those have now been worked out.

interacting_pairs(model, 2.0, :nearest) will match the first nn it scans (if AB = AC, only AB will be chosen) for a given radius
interacting_pairs(model, 2.0, :scheduler) will match nns based on the chosen scheduler (previous implementation, will be useful for interaction hierarchy in multi agent methods for example)
interacting_pairs(model, 2.0, :all) all pairs within a given radius.

[Datseris]

@lhupe I wanted to release version 3.0, but I wanted to wait a bit until there is some progress on your system. I want to test continuous space with more complicated systems which we don't have yet. Do you think you can just open a PR with what you got already and we can help you finish everything?

[lhupe]

My system is more or less finished. However, since the model is based on unpublished work by other people, I will ask for their permission before opening a PR.

[Datseris]

You can also maybe simplify it further so that it isn't connected with actual research?