JuliaDynamics / Agents.jl

Agent-based modeling framework in Julia
https://juliadynamics.github.io/Agents.jl/stable/
MIT License
725 stars 117 forks source link

Bug in periodic continuous spaces due to... floating point error? #853

Closed mastrof closed 1 year ago

mastrof commented 1 year ago

move_agent! and walk! throw position is outside of space extent error, apparently due to floating point approximation

Minimal Working Example

using Agents
pos = (0.4936805014512577, 0.007750092779602313)
vel = (0.6319498548742291, -0.7750092779602327)
dt = 0.01
space = ContinuousSpace((1,1); periodic=true)
model = ABM(ContinuousAgent{2,Float64}, space)
add_agent!(pos, model, vel)
move_agent!(model[1], model, dt)

This particular MWE is on main, but I originally encountered the error on a code running 5.14.0.

What happens is the following (let agent = model[1]): the direction of motion is evaluated as

julia> direction = agent.vel .* dt
2-element SVector{2, Float64} with indices SOneTo(2):
  0.006319498548742291
 -0.007750092779602327

then the "raw" target position would be

julia> target_raw = agent.pos .+ direction
2-element SVector{2, Float64} with indices SOneTo(2):
  0.5
 -1.3877787807814457e-17

which, when passed to normalize_position, becomes

julia> target = Agents.normalize_position(agent.pos .+ direction, model)
2-element SVector{2, Float64} with indices SOneTo(2):
 0.5
 1.0

which errors since 1.0 is not a valid position

mastrof commented 1 year ago

I guess this is a very fringe case if it never came up but I have a set of fairly large simulations that simply cannot be run due to this approximation error.

Tortar commented 1 year ago

Good catch.

What about changing

all(i -> 0 <= pos[i] < space_size[i], 1:D) || error("position is outside space extent!")

in

all(i -> 0 <= pos[i] <= space_size[i], 1:D) || error("position is outside space extent!")

does this have any unwanted side effect?

edit: yes, move_agent! errors since the pos2cell conversion goes outside the grid extent

Tortar commented 1 year ago

But changing also

pos2cell(pos::ValidPos, model::ABM) = Tuple(floor.(Int, pos./abmspace(model).spacing) .+ 1)

to

pos2cell(pos::ValidPos, model::ABM) = Tuple(ceil.(Int, pos./abmspace(model).spacing))

it seems to work fine.

The problem is that this is slightly breaking since for an integer position inside the grid, this will change the result of the position of the agent in the grid

Tortar commented 1 year ago

Actually it is not breaking since the grid is an internal field only used to improve efficiency of neighbors searches. So, is this fine?

mastrof commented 1 year ago

I tried to implement your approach but it breaks a bunch of things, throwing errors of the form

BoundsError: attempt to access 20×20 Matrix{Vector{Int64}} at index [0, 0]
mastrof commented 1 year ago

I guess the problem is because 0 is a valid position in ContinuousSpace. So instead of bumping it up to index 1, ceil(0) returns 0, which is not good for a GridSpace.

mastrof commented 1 year ago

In principle it could be solved by always having a call like normalize_position(pos, abmspace(model).grid) before doing anything with the position returned by pos2cell so it is brought back to a valid grid position, but here I think it gets very tricky because we still want invalid position to throw errors, not being re-mapped into a valid one. E.g. with your correction, you cannot even create an agent at (0,0) because the code cannot place it on the grid; so you would have to add a normalize_position already at creation, but I'm not sure if that would also allow some invalid positions to slip through.

Tortar commented 1 year ago

you are right, then I guess we should just handle the special case with a branch which checks if the value is 0 and makes it 1, hopefully branch prediction should anyway make the code efficient. Otherwise we could clamp it before doing the ceiling

Tortar commented 1 year ago

an even simpler way is to use max between 1 and the ceiling