Implementing the Gillespie algorithm

[MA-Ramirez]

Implementing the Gillespie Algorithm in the agents.jl framework is more complex than I thought. I realised this by coding a new example, the Spatial Rock-Paper-Scissors game. See pull request JuliaDynamics/AgentsExampleZoo.jl#17.

In this issue, I outline the algorithm, the attempted implementation, and the problem. Hopefully, someone comes up with an idea to implement the algorithm, given that it is a common stochastic algorithm used for modelling.

• Gillespie algortihm - Direct Method

1. Initialize the time $t=t_0$ and the system's state $x=x_0$
2. With the system in state $x$ at time $t$, evaluate all the propensities $a_i(x)$ and their sum $a_0(x)$
3. Generate a random value for $\tau$ using an exponential distribution with parameter $a_0(x)$
4. Generate a random value to draw action $i$ in a discrete distribution for which $P(X=i)= a_i/a_0$
5. Update $t \leftarrow t+ \tau$ and $x$ according to $i$
6. Save $(x,t)$ and return to step 1, or else end the simulation

• Attempted implementation

I will focus on steps 2, 3, 4, which are the complex implementation steps.

To take care of steps 2 and 4, I included the following in the agent_step! function:

     #propensities
     a1 = model.selection_rate * nagents(model)
     a2 = model.reproduction_rate * nagents(model)
     a3 = model.exchange_rate * nagents(model)

     #total propensities
     a0 = a1 + a2 + a3

     p = rand(model.rng)
     if 0 <= p < a1/a0
         selection_RPS!(strategy, model)
     elseif a1/a0 <= p < (a1+a2)/a0
         reproduce!(strategy, model)
     elseif (a1+a2)/a0 <= p < 1
         swap!(strategy, model)
     end

For step 3, I define steps as a model property, which is initialized as steps=0. For each model step, the steps increase by one, as follows:

function model_step!(model)
    model.steps += 1
end

The agents have a property called action_step. For each agent, the action step is determined by rand(Exponential(1/a0)). Therefore, when model.steps matches agent.action_step, a given action will occur according to step 4.

• Problem

The values for rand(Exponential(1/a0)) could be of the order of $10^{-5}$ (or any other order less than 1). Therefore the counter model.steps doesn't work to allow the occurrence of actions every time when $t \leftarrow t+ \tau$.

I hope this explanation is helpful.

Cheers, Alejandra R.

[jacobusmmsmit]

If I understand correctly, the actions take place as a continuous time Poisson point process, whereas the steps are incremented by 1 each time. As you identified you would be skipping over a lot of events if you incremented by 1 and events could happen in rapid succession.

It sounds instead like you need an event queue that stores when the next action takes place and by which agent, and have each step fast-forward the model up to this step, tracking the "time" in model.time.

Something like defining a struct to store the events, and using a PriorityQueue from DataStructures.jl:

@enum EventType selection reproduction exchange

struct Event
    what::EventType
    who::SRPSAgent
end

# model pseudocode for initialisation:
model.properties.event_queue = PriorityQueue{Event, Float64}()

then the model step executes the events, and the agent step enqueues the events:

function model_step!(model)
    # 1. Pop and execute event queue
    event, when = dequeue_pair!(model.event_queue)
    model.time = when
    execute_event(event)
    # 2. Do other stuff if you need to
    other_stuff()
end

function agent_step!(model, agent)
    when = rand(abmrng(model), Exponential(1/a0))
    event = Event(whatever_type, agent)
    enqueue!(model.event_queue, event => when)
end

[Datseris]

I have written (on paper) an Agents.jl implementation strategy for Gillespie, because we also need it for a project with @franfram . Currently on the train with bad internet however. Once stable, I'll adjust my implementation to see how it fits with the discussion here and post it here.

[Datseris]

btw welcome @MA-Ramirez thanks for the write up!

[Tortar]

I think that it should be possible to create a new Scheduler which behaves as a PriorityQueue as goodly described by @jacobusmmsmit. With that we can probably generalize this type of procedures

[Datseris]

I don't think this is the most useful way forwards. How would you define the "model step"? How many actions, and/or schedules, should occur within "one step"? I would love to see a design description for this new scheduler because I have trouble seeing it work in practice.

My idea for how to solve this is to create a completely new model type. This is very possible now with (1) the AgentBasedModel API and (2) the change in #889 . Internally the new model has a priority que and agent actions may put more actions into the que. The step! function on this model dispaches either on a function that terminates the evolution, or on a real number, which means to evolve the model for this amount of time and do as many actions as the que will have up to this point. I will be working on this the coming Friday in the meantime @Tortar hasn't shown a concrete design for his proposal.

[jacobusmmsmit]

Perhaps the discrete event simulation (DES) gang could benefit from this too. I for one have written a few DESs and the documentation for these packages in Julia aren't as complete as Agents.jl so the startup cost is too high and I end up writing everything by hand.

Heck, JuliaDynamics itself even has DiscreteEvents.jl (pre-1.0, seems dormant). Given that we are stepping towards API unification with SciML's ecosystem, it would make sense to consider how we could make it easier to perform DES-style simulations nicely either within or combining with Agents.jl.

[Datseris]

https://github.com/JuliaDynamics/ConcurrentSim.jl is also DES and I think more under development.

[Tortar]

Thought about this a bit. I think that indeed a new model could be a good way to implement this feature.

However I have a consideration to share: this model doesn't exclude the others in the sense that we could have this new model either with inside either a vector or a dict for the agents. But this means (correct me if I'm wrong) that we will actually have a proliferation of model types since some can be mixed with each other, so that we actually need an "UnremovableEventBasedModel" and a "StandardEventBasedModel" (what about when we will have a "MultiAgentModel" optimized for multi-agents models?). So I'm thinking if this kind of separation in different types for the model is optimal or not, maybe it would be better to have only one unique exported function for the model (e.g. ABM(....) with keywords to choose the actual internal representation) and then manage things internally in the simplest way, which can be with either new types or dispatch on the parameters of the struct? What do you think?

[Datseris]

I disagree with ABM(...) doing different things based on keywords. Functions and names should convey the purpose of the thing and ultra generic names should be avoid in favor of specific ones. Leads to more clarity and less confusion and less likely to do a mistake. Not only that, but "dispatch on keywords" is not a core language feature, in contrast to multiple dispatch, which is native in the language.

The problem you raise can be addressed by allowing users to configure the container type and having internal functions for adding and removing agents from the container type, irrespectively of the over-arching model type. In fact, we already have these functions. They are called add_to_model!. They are each 2 lines of code. And we would need two more lines of code for each version for the new Gillespie model. These extra 4 lines of code are not enough (in my eyes) to justify the ambiguity that will come with simply naming everything ABM.

[Tortar]

The configuration of the container type by the users is indeed one of the solutions I thought about. I have also an idea about how to allow removals in a UnremovableABM without impacting performance of the currently supported methods so that both containers could share the same set of operations.

[Tortar]

Another thing worth considering is if we can add a new component to SingleContainerABM to trigger the dispatch on step! (probably a new scheduler is not enough, but some other sort of component can be). We can maybe avoid the new model type this way and just add the component to the already existing model types

[Datseris]

I don't see a reason to avoid the new type. It will only increase clarity to have different types with different rules. Gillespie models are stepped over real time. Standard ABM is stepped over discrete time. Why fixate so much in keeping one type? It will anyways be impossible once we figure out a solution for multi agent models.