Simple discrete event simulation in Julia

I am porting a process-based discrete simulation DSL from Racket into Julia and am beginning by developing a very simple (almost trivial) sample to make sure I understand using tasks to represent simulation processes and can schedule and execute them properly.

The code below does work and is about ten times faster than the corresponding Racket code. But, I was wondering if others may have suggestions to improve it - particularly in the tasking area (e.g., alternatives to yieldto). I am not real happy with the yieldto I had to put in the simulation processes (generator and customer) but they will be hidden inside a @process macro eventually. Stylistic comments, etc. are more than welcome.

I also apologize in advance if posting code of this nature or length is frowned upon - this is my first post here. If this isn’t the desired way to do it, please let me know the preferred method.

Thanks in advance
Doug Williams

# Simplified Simulation System

using Printf
using Random
import Base.schedule

"Specify whether the simulation execution should be traced."
const trace = false

# Event Definition and Scheduling

"""
An event represents the application of procedure to its arguments at some
(simulated) future time.
"""
struct Event
    time::Float64
    task::Task
end

isgreater(x::Event, y::Event) = x.time > y.time

"""
A simulation represents the state of an execution. Specifically, it contains
the current time, current event, future event list, and control task for the
simulation.
"""
mutable struct Simulation
    current_time::Float64
    current_event::Union{Event, Nothing}
    future_event_list::Array{Event, 1}
    control_task::Task
    Simulation() = new(0.0, nothing, [], current_task())
end

"The current active simulation."
const current_simulation = Simulation()
current_time() = current_simulation.current_time
current_event() = current_simulation.current_event
control_task() = current_simulation.control_task

"""
    event_schedule(ev::Event, ev_list::Array{Event, 1})
    
Adds an event the the given event list. This routine keeps the event list in
sorted (ascending) order.
"""
function event_schedule(ev::Event, ev_list::Array{Event, 1})
    index = searchsortedfirst(ev_list, ev, lt=isgreater)
    insert!(ev_list, index, ev)
end

"""
    schedule(ev::Event)
    schedule(sim::Simulation, ev::Event)

Schedules the event on the future event list for the specified simulation,
which defaults to the current simulation.
"""
schedule(ev::Event) = schedule(current_simulation, ev)
function schedule(sim::Simulation, ev::Event)
    event_schedule(ev, sim.future_event_list)
end

# Simulation control routines

"""
    work(delay::Real)

Simulate the delay while work is being done.  Add an event to return to this task
in the future to the event list.
"""
function work(delay::Real)
    schedule(Event(current_time() + delay, current_task()))
    yieldto(current_simulation.control_task)
end

"""
    start_simulation()

This is the main simulation loop.
"""
function start_simulation()
    while !isempty(current_simulation.future_event_list)
        current_simulation.current_event = pop!(current_simulation.future_event_list)
        current_simulation.current_time = current_simulation.current_event.time
        if trace
            @printf("%9.3f: %s\n", current_time(),
                    current_simulation.current_event.task)
        end
        yieldto(current_simulation.current_event.task)
        current_simulation.current_event = nothing
    end
end

# Random Distributions

"Returns a random float from a uniform distribution between a and b."
function random_uniform(a::T, b::T)::T where T <: AbstractFloat
    a + rand(T) * (b - a)
end

"Returns a random float from an exponential distribution with mean mu."
function random_exponential(mu::T)::T where T <: AbstractFloat
    -mu * log(rand(T))
end

# Example simulation model

# "Process to generate n customers arriving into the system."
# @process function generator(n::integer)
#     for i = 1:n
#         work(random_exponential(4.0))
#         @schedule now customer(i)
#     end
# end

"Process to generate n customers arriving into the system."
function generator(n::Integer)
    @task begin
        try
            for i = 1:n
                work(random_exponential(4.0))
                schedule(Event(current_time(), customer(i)))
            end
        finally
            yieldto(current_simulation.control_task)
        end
    end
 end

# "The ith customer into the system."
# @process function customer(i::Integer)
#     println("$(current_time()): customer $i enters")
#     work(random_uniform(2.0, 10.0))
#     println("$(current_time()): customer $i leaves")
# end

"The ith customer into the system."
function customer(i::Integer)
    @task begin
        try
            println("$(current_time()): customer $i enters")
            work(random_uniform(2.0, 10.0))
            println("$(current_time()): customer $i leaves")
        finally
            yieldto(current_simulation.control_task)
        end
    end
end

# "Run the simulation for n customers."
# function run_simulation(n::Integer)
#     @with_new_simulation begin
#         @schedule at 0.0 generator(n)
#         start_simulation()
#     end
# end

# "Run the simulation for n customers."
function run_simulation(n::Integer)
    schedule(Event(0.0, generator(n)))
    start_simulation()
end

@time run_simulation(10)
1 Like

You probably don’t need random_uniform; rand(1.5:2.5) samples from a uniform distribution. This is available immediately, without Distributions.jl. You might consider using Distributions instead of random_exponential, especially if want to use more distributions.

It looks to me like event_schedule orders the array by decreasing time, not increasing. To get increasing time, you can define a method for Base.isless and Base.isequal. Then you don’t need to pass the keyword arg lt=....

Most Julia style guides suggest using an explicit return statement. I like to follow this rule. Sometimes it’s not clear if a big statement at the end of a function is only there for its side effect.

I think this samples from the elements of the range 1.5:2.5, i.e. you get 1.5 half the time and 2.5 half the time. I think random_uniform is needed here (unless you add Distributions as a dependency).

2 Likes