GPU Extensions to Agents.jl based on a different modelling approach

[gatocor]

I developed a package for agent based modelling in continuum spaces of 1D, 2D and 3D.

dsb-lab/AgentBasedModels.jl

The features that it had over the current implementation of Agents.jl are:

1. GPU compatible.
2. Implemented continuum model integration methods as Euler, Heun or RungeKutta4.
3. The possibility to define a continuum field coupled with the agents.
4. Optimization methods specific for non-differentiable problems as agent-models.

Although the package is more specific purposed than Agents.jl. Talking with @Datseris, we agreed that it include interesting features and we could try to make it an extension or even a merge into Agents.jl.

On the compatibility with GPU

Since there are quite a lot of different topics, we can start talking about the GPU compatibility of Agents.

So far, Agents is arranged as:

• Array of mutable structures (agents).
• Model with the global parameters of the model

In order for Agents.jl to be compatible with CUDA.jl, it is required to transform this vision as:

• Each parameter of an agent to become a CUDA.Array by itself.
• Model may contain both global and local parameters.

This transformation will not need mutable structures that are incompatible with CUDA.jl and allow GPU kernels to work out.

On the way we define code in the other package

For allowing models to work in both CPU and GPU, my package uses macro programming to vectorize the code in the appropriate shape compatible to CPU and GPU. It will be interesting to discuss how we could make compatible Agents.jl way of defining stepping functions to allow for this flexibility.

[Datseris]

(let's keep this issue dedicated to the GPU discussion. Other functionality of the package seems orthogonal to the GPU and could be discussed in another issue. I've changed the title to reflect this)

---

Can you please paste here the code you've showed me in the Jupyter notebook during our call, that had a trivial model that (1) updated an agent property (2) moved an agent and (3) did some operation based on nearby agents of an agent?

[gatocor]

Here it is a small model with a random walker that basically counts how many neighbors it has at distance below 0.1.

Model definition:

model = Agent(2, #Dimensions
    
    #Local parameters
    localIntInteraction = [:nn], #Parameter that is reset to 0 at every step
    
    #Global parameters
    globalFloat = [:D], #Global parameter of the diffusion process
    
    #Agent rules
    updateInteraction = quote
       
        if euclideanDistance(x.i, x.j, y.i, y.j) < 0.1
            
            nn += 1            
            
        end
        
    end,
        
    updateLocal = quote
       
        x += Normal(0,D)*dt
        y += Normal(0,D)*dt
        
    end,
    
    neighbors = :CellLinked, #Method for computing neighbors
    platform = :CPU #Platform in which to evaluate the neighbors
);

When created, the model creates using macros adapting the code to the platform and saving them compiled in the model object. You can see how they look like in model.declaredUpdatesCode

> prettify(model.declaredUpdatesCode[:UpdateLocal_])

:(
(t, dt, N, NMedium, nMax_, id, idMax_, 
simBox, NAdd_, NRemove_, NSurvive_, 
flagSurvive_, holeFromRemoveAt_, repositionAgentInPos_, 
skin, dtNeighborRecompute, nMaxNeighbors, cellEdge, 
flagRecomputeNeighbors_, flagNeighbors_, neighborN_, 
neighborList_, neighborTimeLastRecompute_, posOld_, 
accumulatedDistance_, nCells_, cellAssignedToAgent_, 
cellNumAgents_, cellCumSum_, x, y, z, 
xNew_, yNew_, zNew_, varAux_, varAuxΔW_, 
liNM_, liM_, liMNew_, lii_, lfNM_, lfM_, lfMNew_, lfi_, 
gfNM_, gfM_, gfMNew_, gfi_, giNM_, giM_, giMNew_, gii_, 
mediumNM_, mediumM_, mediumMNew_) #Arguments that are provided to all functions, most are empty in this simple case
->
@inbounds(Threads.@threads(for i1_ = 1:1:N[1]
                  xNew_[i1_] += AgentBasedModels.rand(AgentBasedModels.Normal(0, gfNM_[1])) * dt[1]
                  yNew_[i1_] += AgentBasedModels.rand(AgentBasedModels.Normal(0, gfNM_[1])) * dt[1]
              end))
)

Creation of a Community of agents:

com = Community(model,
#Setting the number of agents
N=[10],
#Parameters required for the neighbor algorithm
simBox=[-10 10;-10 10.],
cellEdge=[2,2.],
#Step size of the time
dt=[0.1],
)

#Initialize the global parameter externaly (can also be done internally)
com.D = 0.1

The community has the information of all the agent local and global parameters in vector format:

> com.x
[0,0,0,0,0,0,0,0,0,0]
> com.nn
[0,0,0,0,0,0,0,0,0,0]

Evolving all the set of rules for some number of steps:

evolve!(com,steps=100)

[AayushSabharwal]

This looks really cool. At a glance, it strikes me as similar in architecture to https://github.com/cesaraustralia/DynamicGrids.jl, where instead of the grid we have (from a user perspective) a vector of agents. I'll have to dive more deeply to see if it really is all that similar, which I'll get around to soon.