Sample Main Loop


Let's walk through a small simulation so that we can see how CompetingClocks could appear in the main loop. This models individuals wandering around on a checkerboard, where no two individuals can occupy the same space. First, let's import libraries.

import Base
using Distributions
using Random
using SparseArrays
using Test
using CompetingClocks

There is a physical state for the model, separate from the state of the bag of clocks and separate from the time. Most of the board is empty, so let's use a spare matrix to represent the locations of individuals.

mutable struct PhysicalState
    board::SparseMatrixCSC{Int64, Int64}

# They can move in any of four directions.
@enum Direction Up Left Down Right
const DirectionDelta = Dict(
    Up => CartesianIndex(-1, 0),
    Left => CartesianIndex(0, -1),
    Down => CartesianIndex(1, 0),
    Right => CartesianIndex(0, 1),

The simulation, itself, carries the state of the clocks in the sampler, as well as the physical state. We'll call it a finite state machine (FSM) because it has the traits of a Moore machine.

mutable struct SimulationFSM{Sampler}

Main Loop

The main loop will ask what event happens next to the state. When that event changes the state, the loop will update the set of possible next events by disabling outdated ones and enabling new ones. The calls to CompetingClocks are:

  • next(sampler, current time, random number generator (RNG))
  • enable!(sampler, event ID, distribution, current time, start time of clock, RNG)
  • disable!(sampler, event ID, current time)

There are a lot of samplers in CompetingClocks to choose from. This example uses CombinedNextReaction algorithm, which has good performance for a variety of distributions. Samplers in CompetingClocks require two type parameters, a key type for clocks and the type used to represent time. In this case, the clock key type fully represents an event, giving the ID of the individual, where they start, and which direction they may move.

const ClockKey = Tuple{Int,CartesianIndex{2},Direction}

function run(event_count)
    Sampler = CombinedNextReaction{ClockKey,Float64}
    physical = PhysicalState(zeros(Int, 10, 10))
    @test showable(MIME("text/plain"), physical)
    sim = SimulationFSM{Sampler}(
    initialize!(sim.physical, 9, sim.rng)
    current_events = allowed_moves(sim.physical)
    for event_id in current_events
        enable!(sim.sampler, event_id, Weibull(1.0), 0.0, 0.0, sim.rng)

    for i in 1:event_count
        (when, what) = next(sim.sampler, sim.when, sim.rng)
        if isfinite(when) && !isnothing(what)
            sim.when = when
            move!(sim.physical, what)
            next_events = allowed_moves(sim.physical)
            for remove_event in setdiff(current_events, next_events)
                disable!(sim.sampler, remove_event, when)
            for add_event in setdiff(next_events, current_events)
                enable!(sim.sampler, add_event, Weibull(1.0), when, when, sim.rng)
            current_events = next_events
            @show (when, what)

For this checkerboard with wandering individuals, the allowed moves are moves by any individual to any free square on the board. The set of allowed moves is precisely the set of enabled clocks, so it stores ClockKeys.

function allowed_moves(physical::PhysicalState)
    allowed = Set{ClockKey}()
    row, col, value = findnz(physical.board)
    for ind_idx in eachindex(value)
        location = CartesianIndex((row[ind_idx], col[ind_idx]))
        for (direction, offset) in DirectionDelta
            if checkbounds(Bool, physical.board, location + offset)
                if physical.board[location + offset] == 0
                    push!(allowed, (value[ind_idx], location, direction))
    return allowed

Changes to the state of the board

We set up the board with an initializer that places individuals at random. We move one individual at a time, when their next event happens.

function initialize!(physical::PhysicalState, individuals::Int, rng)
    physical.board .= 0
    locations = zeros(CartesianIndex{2}, individuals)
    for ind_idx in 1:individuals
        loc = rand(rng, CartesianIndices(physical.board))
        while physical.board[loc] != 0
            loc = rand(rng, CartesianIndices(physical.board))
        locations[ind_idx] = loc
        physical.board[loc] = ind_idx

function move!(physical::PhysicalState, event_id)
    (individual, previous_location, direction) = event_id
    next_location = previous_location + DirectionDelta[direction]
    # This sets the previous board value to zero.
    SparseArrays.dropstored!(physical.board, previous_location.I...)
    physical.board[next_location] = individual

How it runs

A run of this simulation produces a sequence of moves, no two happening at the same time because this is a continuous-time simulation.

(when, what) = (0.07319364933011555, (8, CartesianIndex(2, 10), Main.var"##381".Down))
(when, what) = (0.10866577949385807, (2, CartesianIndex(5, 10), Main.var"##381".Down))
(when, what) = (0.16874877793434204, (6, CartesianIndex(8, 8), Main.var"##381".Down))
(when, what) = (0.18621506968095888, (6, CartesianIndex(9, 8), Main.var"##381".Down))
(when, what) = (0.20347320994043128, (3, CartesianIndex(1, 8), Main.var"##381".Down))
(when, what) = (0.20635380502472672, (6, CartesianIndex(10, 8), Main.var"##381".Right))
(when, what) = (0.2118081725791219, (5, CartesianIndex(4, 1), Main.var"##381".Right))
(when, what) = (0.2634608206635203, (1, CartesianIndex(7, 10), Main.var"##381".Left))
(when, what) = (0.26949661650175905, (3, CartesianIndex(2, 8), Main.var"##381".Right))
(when, what) = (0.3035113012060176, (2, CartesianIndex(6, 10), Main.var"##381".Left))

This page was generated using Literate.jl.