Modelling collisions of agents with the environment in discrete time

Specific Domains Modelling & Simulations
1

Hi everyone,

I am new to agent-based modelling and have recently started using Agents.jl for a research project. Thanks a lot for this great package!

The model should simulate how (for example) coconuts disperse through the pacific ocean with the help of ocean currents.

One of the goals is to track whether / on which islands the coconuts arrive. I thought I could achieve this by adding an island_matrix to the model properties, and removing an agent when its location is marked as an island.

The problem I have run into here is that – depending on the current’s strength – agents “hop over” islands, leading to the contact with the island never being recorded and the agent remaining in the simulation.

For an illustration, please see this MWE:

using Agents, Random, CairoMakie

@agent struct Cocos(ContinuousAgent{2,Float64})
    lifetime::Int
    pos_type::Int
end

function initialize(;
    n_nuts=14,
    lifetime=40,
    extent=(150, 150),
    seed=1,
)
    space2d = ContinuousSpace(extent; spacing=1.0, update_vel!)
    rng = Xoshiro(seed)

    xdim, ydim = extent

    island_matrix = zeros(extent)
    island_matrix[:, 100] .= 1

    properties = (
        current=[(0, 1) for _ in 1:xdim, _ in 1:ydim],
        island_matrix=island_matrix,
    )

    model = StandardABM(Cocos, space2d; properties, rng, agent_step!, model_step!, scheduler=Schedulers.Randomly())

    for i in 1:n_nuts
        vel = (0, 0)
        pos = (10i, ydim / 2)
        pos_type = getindex_grid(model.island_matrix, pos)
        add_agent!(pos, model, vel, lifetime, pos_type)
    end

    return model
end

function agent_step!(agent, model)
    if expired(agent) || onisland(agent)
        remove_agent!(agent, model)
    else
        # illustrate the effect of a strong current
        # in the real model, `dt` will be 1.0 and
        # fast movement would result from values in `model.current`
        dt = 10.0
        move_agent!(agent, model, dt)
    end
    return nothing
end

function model_step!(model)
    # more stuff here in the real model ...
    update_agents!(model)
end

function update_vel!(agent, model)
    agent.vel = getindex_grid(model.current, agent.pos)
end

function update_agents!(model)
    for agent in allagents(model)
        agent.lifetime -= 1
        maybe_update_pos_type!(agent, model)
    end
    return nothing
end

function maybe_update_pos_type!(agent, model)
    if onisland(agent, model)
        agent.pos_type = getindex_grid(model.island_matrix, agent.pos)
    end
    return nothing
end

onisland(agent::Cocos, model) = getindex_grid(model.island_matrix, agent.pos) > 0

onisland(agent::Cocos) = agent.pos_type > 0

expired(agent::Cocos) = iszero(agent.lifetime)

truncate(x; at=1) = x < at ? 1 : x

coords(v) = truncate.(round.(Int, v, RoundUp))

getindex_grid(A, pos) = A[coords(pos)...]

model = initialize()

const cocos_polygon = Makie.Polygon(Point2f[(-1, -1), (2, 0), (-1, 1)])

function cocos_marker(c::Cocos)
    φ = atan(c.vel[2], c.vel[1])
    rotate_polygon(cocos_polygon, φ)
end

islandcolor(model)::Matrix{Int} = model.island_matrix .> 0

heatkwargs = (colormap = [:white, :black], colorrange = (0, 1))

plotkwargs = (;
    agent_marker = cocos_marker,
    agent_color = :orange,
    agentplotkwargs = (strokewidth = 1.0, strokecolor = :black),
    heatarray = islandcolor,
    heatkwargs = heatkwargs,
)

abmvideo(
    "floating.mp4", model;
    framerate=10, frames=40,
    title="Coconuts",
    plotkwargs...
)

I’m sure that there are other oddities in my code; for example, it felt a bit strange to index the ocean current force field having to first convert agent.pos (a Float64) to Int, and then ensuring that both coordinates are > 0. Is this how it’s supposed to be done?

Thank you for your help in advance! :slight_smile:

1 Like
2

I’ve found a (somewhat imprecise) solution that is acceptable for my use case. In short, this solution tracks the last position of each agent, and at each step checks whether there is an island in the space spanned by the previous and current positions.

This solution is great when an agent moves only along a single direction, but gets less precise the more diagonal the movement. Since the movement in my application will be quite slow and the spatial resolution fairly high, I think that these imprecisions won’t matter.

using Agents, Random, CairoMakie

@agent struct Cocos(ContinuousAgent{2,Float64})
    lifetime::Int
    pos_type::Int
    last_pos::SVector{2, Float64}
end

function initialize(;
    n_nuts=14,
    lifetime=40,
    extent=(150, 150),
    seed=1,
)
    space2d = ContinuousSpace(extent; spacing=1.0, update_vel!)
    rng = Xoshiro(seed)

    xdim, ydim = extent

    island_matrix = zeros(Int, extent)
    island_matrix[:, 100] .= 1

    properties = (
        current=[(0, 1) for _ in 1:xdim, _ in 1:ydim],
        island_matrix=island_matrix,
    )

    model = StandardABM(Cocos, space2d; properties, rng, agent_step!, model_step!, scheduler=Schedulers.Randomly())

    for i in 1:n_nuts
        vel = (0, 0)
        pos = (10i, ydim / 2)
        pos_type = getindex_grid(model.island_matrix, pos)
        last_pos = pos
        add_agent!(pos, model, vel, lifetime, pos_type, last_pos)
    end

    return model
end

function agent_step!(agent, model)
    if expired(agent) || onisland(agent)
        remove_agent!(agent, model)
    else
        agent.last_pos = agent.pos
        move_agent!(agent, model, 5.0)
    end
    return nothing
end

function model_step!(model)
    update_agents!(model)
end

function update_vel!(agent, model)
    agent.vel = getindex_grid(model.current, agent.pos)
end

function update_agents!(model)
    for agent in allagents(model)
        agent.lifetime -= 1
        maybe_update_pos_type!(agent, model)
    end
    return nothing
end

function maybe_update_pos_type!(agent, model)
    if onisland(agent, model)
        agent.pos_type = in_path(model.island_matrix, agent.last_pos, agent.pos)
    end
    return nothing
end

function onisland(agent, model)
    island_in_path = in_path(model.island_matrix, agent.last_pos, agent.pos)
    return !isempty(island_in_path)
end

onisland(agent::Cocos) = agent.pos_type > 0

expired(agent::Cocos) = iszero(agent.lifetime)

truncate(x; at=1) = x < at ? 1 : x

coords(v) = truncate.(round.(Int, v, RoundUp))

getindex_grid(A, pos) = A[coords(pos)...]

sorted_bounds(a, b) = a < b ? (a, b) : (b, a)

function in_path(env, a, b)
    xa, ya = coords(a)
    xb, yb = coords(b)

    # Exit if coordinates match after rounding
    if xa == xb && ya == yb
        return 0
    end
    
    x_min, x_max = sorted_bounds(xa, xb)
    y_min, y_max = sorted_bounds(ya, yb)
    
    # Search for island but error if more than one island in path
    island_tile = 0
    @inbounds for x in x_min:x_max
        for y in y_min:y_max
            current_tile = env[x, y]
            if current_tile != 0
                if island_tile == 0
                    island_tile = current_tile
                elseif island_tile != current_tile
                    error("Two or more islands in path: [$current_tile, $island_tile]")
                end
            end
        end
    end
    
    return island_tile
end

model = initialize()

const cocos_polygon = Makie.Polygon(Point2f[(-1, -1), (2, 0), (-1, 1)])

function cocos_marker(c::Cocos)
    φ = atan(c.vel[2], c.vel[1])
    rotate_polygon(cocos_polygon, φ)
end

islandcolor(model)::Matrix{Int} = model.island_matrix .> 0

heatkwargs = (colormap = [:white, :black], colorrange = (0, 1))

plotkwargs = (;
    agent_marker = cocos_marker,
    agent_color = :orange,
    agentplotkwargs = (strokewidth = 1.0, strokecolor = :black),
    heatarray = islandcolor,
    heatkwargs = heatkwargs,
)

abmvideo(
    "floating.mp4", model;
    framerate=20, frames=15,
    title="Coconuts",
    plotkwargs...
)