I am working to improve the performance of an agent-based model. The model is of two bacterial species whose growth and death probabilities are determined by the amount of a nutrient available locally. I’m aware of Agents.jl but I haven’t had a chance to learn it so this model is basically just written in base Julia.
The domain of the model is two parallel k x k grids. One grid contains bacteria and the other contains a nutrient. In the bacteria grid, spaces are sampled until a certain number of bacteria have been selected per time step. When selected, a bacteria can either reproduce, die, or do nothing, with probabilities determined by the amount of nutrient in the corresponding space in the second grid. The amount of nutrient in each space of that grid is determined by solving an ode in that space over one time step of the agent-based model. So the nutrient grid determines a kxk ode system that can be either coupled or decoupled depending on parameter values.
I also have a spatially homogenous version of the model where bacteria can place offspring anywhere in the spatial domain. In this case, there’s a global amount of nutrient that determines growth/death probabilities and only one ode that needs to be solved each time step. For a 40x40 grid size, the homogenous version takes ~3 seconds on my machine and the heterogenous version 3-4 minutes.
So the problem in the heterogeneous version is clearly in solving the ode system each time step. There are also some functions that do things like check if a bacteria is on a boundary and add up all the nutrient among neighboring spaces that add some cost. What I’m currently doing for the ode system is looping over the entire grid, adding up nutrient, and defining and solving an ode problem. I think this is where there could be room for improvement but I’m not sure what would work.
Would it be better to post a full MWE (~280 lines) or just sections of my code that I think are causing the issue?
Why use an ODE here? Why not a simple discrete time-stepping rule? ODE solvers tend to do a lot of work to get answers accurate to 8 significant figures or the like. Is that necessary in this simulation?
Can you explain a discrete-time stepping rule a little more? In the past, I had just been using a hand-coded RK4 to advance the odes, but this got to be an issue when I put in realistic parameter values. Some parameters are very large so I switched to the ode solver to handle stiffness issues
Well, can you give the definition of the ODE? then we could maybe deconstruct the model a bit, see what things are really needed and what might be simply approximated.
I was also wondering why you need an ODE. A discrete grid with rules updating state over many timesteps is already like a fixed timestep ODE (or lots of them). So youre nesting ODEs… You can instead reduce the time step if you need more accuracy.
There is also DynamicGrids.jl where you can run something like this parallel or on GPUs, with pretty clean code. A few orders of magnitude faster than Agents.jl for things they cab both do, but is designed for working with populations and quantities rather that individually mobile/tracked individuals (Until I write a GPU compatible mobile agent layer).
where lambda, mu, and eta are rate parameters, C is 0 if the space is empty and 1 if a bacteria is there. m is the number of neighboring spaces there are (m = 8 for most spaces, fewer on a boundary space), g is a diffusion rate between neighbors, and X is the sum of all the nutrients in the neighbor spaces.
I think when I was using RK4 to that was basically a discrete-time update and I was choosing a step size as large as possible to get reasonable speed. Messing with the step-size settings in DifferentialEquations.jl is something I still need to try. But I don’t know how much of my issue is just do to having a ton of loops vs time spent updating the nutrient grid (ODE or no…)
To me it seems like you don’t need DifferentialEquations.jl at all here, its just complicating the problem. You just need a simple fixed time-step simulation?
If bacteria is presence/absence only, you can do this with DynamicGrids.jl. It will be (many) orders of magnitude faster than this. That should probably run in a tiny fraction of a second, not minutes. You propably also want to wrap the boundaries as a torus instead handling complex boundary conditions like that. DynamicGrids.jl handles that for you.
Yes, DynamicGrids.jl is what I’m looking for. Basically I just need to do two grids and have rules for births, deaths, and updating the second grid? And rules can be written such that they take inputs from both grids?
Pretty much. You can take input from as many grids as you like. I imagine something like this is what you want. It’s not actually checked, and I’ve swapped C to b to represent values from the bacteria grid, and x to n to represent values from the nutrient grid :
nutrient_rule = let lambda=lambda, mu=mu, eta=eta
Neighbors{Tuple{:nutrient,:bacteria},:nutrient}}(Moore(1)) do data, hood, (n, b), I
X = sum(hood) # amount of nutrient in neighborhood
n + lambda - mu*n - eta*b*x - g*m*n + g*X,
end
end
bacteria_rule = let somparam=someparam
Neighbors{Tuple{:bacteria,:nutrient},:bacteria}}(Moore(1)) do data, hood, (b, n), I
sum(hood) # number of bacteria in neighborhood
# some birth/death rule
end
end
init = (bacteria=rand(Bool, 40, 40), nutrient=rand(40, 40))
# You could also use a visual output here.
# E.g. InteractOutput will let you control the parameters
output = ArrayOutput(
init=init,
tspan=1:100 # Or use Unitful or Dates
)
sim!(output, bacteria_rule, nutrient_rule; boundary=Wrap())
The let blocks are for performance, to avoid global variables.
The main examples in the docs are all for single grid rules, which are easier to understand, but don’t really demonstrate the more complicated models. I should add a model zoo for multi-grid rules. If you end up with something working we could add it as an example.
How are the rules applied each time-step? In my code in each time step, I sample the bacteria grid until it encounters N individuals, where N was the total population at the beginning of the time step.
If you happened to have any examples of multi-grid rules that would be helpful. I also have rules about where offspring bacteria are placed that right now are done with a lot of confusing conditional logic, so I’m wondering how complicated the rules in DynamicGrids.jl can/should be?
Rules are applied to every cell in the grid, for every rule in order. Then the timestep is incremented, and it all runs again. There is array swapping and boundary handling happening in the background, but you can mostly ignore that.
I sample the bacteria grid until it encounters N individuals, where N was the total population at the beginning of the time step.
I’m not sure I understand what you mean here. What does “I sampled” mean, and what is “it” in “until it”? But if I do - you mean you are running your rule for the first N bacteria only? Then why not just run for all of the bacteria? This should be very fast on a 40x40 grid.
And you can’t actually run a subset of cells in DynamicGrids.jl, it breaks an important condition: rules running in the same timestep don’t know anything about each other, and there is no known order of application accross the grid. We are creating the illusion that all rule applications happen at exactly the same time, and state changes are only visible in the next time step. And they might actually all run at the same time on a GPU.
I’m wondering how complicated the rules in DynamicGrids.jl can/should be?
DynamicGrids.jl models can get complicated. I’ve used e.g. ~10 400*400 grids, ~15 rules, although that’s getting a bit much. But part of the goal of DynamicGrids.jl is making code for something so complicated easy to understand and modular.
If you mean how complicated should individual rules be, It really depends on what you have to do. This is the most complicated I know of, and that’s mostly setup code. But normally individual rules are under 10 lines.
What I mean is, let’s say there are N bacteria at the start of a time step. Over a time step, I randomly sample cells in the grid. If a cell is empty then nothing happens and pick another cell, and if a cell is occupied then rules are applied and a counter is incremented. Then cells are randomly sampled until the counter gets to N.
I *believe* the reason for this is because the number of occupied cells changes during the course of a time step and this keeps each time step consistent in some way, but I might not know a rigorous reason for it.
But I’m just trying to understand if I’m doing the same thing as DynamicGrids.jl and if I am then I should be good
I think randomly sampling until the counter hits N is a bad method. You’ll likely resample some cells, and miss others, and it’ll probably take longer than just going through all the cells.
You’re doing a very different thing to DynamicGrids.jl, but also to most simulations. I would like to understand the reason a little more. Normally population change is meant to affect the dynamics… otherwise why have a birth/death rule at all.
Do you want a moving population that doesn’t change in number? If so, I would just write the simulation like that. Moving bacteria that don’t die. You can either use a grid of Int so each cell can hold multiple bacteria, then just randomly move them around, or write something more complicated where they only move to unnocupied cells. Use a SetNeighborhoodRule so you can move them and remove the source cell manually, or e.g. not move them if there is no unnocupied cell. Its the opposite of a Neighborhood rule you can write to the neighborhood instead of reading from it to change the current cell.