Discrete event simulation framework for discrete times in Julia

Question

I would like to receive constructive feedback on my simulation framework (this is for my thesis work). I have no formal background in programming, and I am wondering if there is a more effective and efficient way to structure my codes for better readability and scalabiilty so that people can easily understand, modify, and build on my codes while achieving high computational efficiency.

Here is a brief description of my simulation model. For each agent I create in my simulation world, they go through a number of event processes annually until they die, reach the maximum age (e.g., 111), or reach the max. time window (e.g., 2050 if an agent is created in 2015 and I set the time window to be 35). I am using a modular apporach, in which each module contains at least one process (so far I have not encoutered more than 1, but I think I would in the future): a process function for that module and input parameters for that process function.

To facilliate that framework, I have created an abstract type Simulation_Module under which event modules (birth and death) are contained.

abstractModule.jl:

abstract type Agent_Module end

abstract type Simulation_Module end

# sub modules
abstract type Birth_Module <: Simulation_Module  end
abstract type Death_Module <: Simulation_Module  end

agent.jl: Here is an agent module

mutable struct Agent{Z<:Bool,Y<:Int64}  <: Agent_Module
    sex::Z
    age::Y
    birth_year::Y
    cal_year::Y
    alive::Z
end

function set_agent!(ag, sex, age, birth_year,cal_year,alive)
    ag.sex = sex;
    ag.age= age;
    ag.birth_year = birth_year;
    ag.cal_year = cal_year;
    ag.alive = alive;
    nothing
end

agent = Agent(false,0,2015,2015,true)

birth.jl: Birth module

mutable struct Birth <: Birth_Module
    parameters::Union{AbstractDict,Nothing}
    process::Function
end

parameter_names = nothing

function birth_process(cal_year_index::Int64,parameters::Union{AbstractDict,Nothing}=nothing)::Bool
    rand(Bernoulli(0.5))
end

birth = Birth(Dict_initializer(parameter_names),birth_process)

death.jl

mutable struct Death <: Death_Module
    parameters::Union{AbstractDict,Nothing}
    process::Function
end

parameter_names = nothing

function death_process(ag::Agent,parameters::Union{AbstractDict,Nothing}=nothing)::Bool
    rand(Bernoulli(0.5))
end

death = Death(Dict_initializer(parameter_names),death_process)

Simulation.jl: It contains a struct for Simulation and a run function for running a simulation.

using Setfield, Distributions, StatsFuns
using TimerOutputs

function Dict_initializer(parameter_names::Union{Nothing,Vector{Symbol}})
    isnothing(parameter_names) ? nothing : Dict(parameter_names .=> missing)
end

# abstract types
include("abstractModule.jl")
# agent
include("agent.jl")
# event proceses
include("birth.jl")
include("death.jl")


mutable struct Simulation <: Simulation_Module
    max_age::Int64
    starting_calendar_year::Int64
    time_horizon::Int64
    n::Int64 # number of agents to create in each year
    Agent::Agent_Module
    Birth::Birth_Module
    Death::Death_Module
    OutcomeMatrix::AbstractDict # needs to be a dictionary
end


function run(simulation::Simulation_Module)

    # Is it better to define the fields of the simulation object
    # separately if I am going to use
    # them in simulation multiple times?
    max_age::Int64 = simulation.max_age
    min_cal_year::Int64 = simulation.starting_calendar_year
    max_cal_year::Int64 = min_cal_year + simulation.time_horizon - 1
    simulation.n = (isnothing(simulation.n) ? 100 : simulation.n)
    n::Int64 = simulation.n

    max_time_horizon::Int64 = simulation.time_horizon
    cal_years::Vector{Int64} = collect(min_cal_year:1:max_cal_year)

    # store events
    # total num of agents
    n_list = zeros(Int64,simulation.time_horizon,2)
    # store events by cal year, age, and sex for each event type
    event_list::Vector{String} = ["death","alive"]
    event_matrix = zeros(Int64,length(cal_years),max_age,2,length(event_list))

    # time the performance
    to = TimerOutput()

    # initiate variables
    tmp_n::Int64 = 0
    tmp_cal_year_index::Int64 = 0

    @timeit to "sleep" sleep(0.02)
    for cal_year in cal_years
        # for each calendar year
        @timeit to "calendar year $cal_year" begin
            tmp_cal_year_index = cal_year - min_cal_year + 1

            for i in 1:n
                # initialization: create a new agent
                set_agent!(simulation.Agent,simulation.Birth.process(tmp_cal_year_index),0,tmp_cal_year_index,tmp_cal_year_index,true)

                # update the total number of agents added to the simulation
                n_list[tmp_cal_year_index,simulation.Agent.sex+1] +=1

                # go through event processes for each agent until
                # death or max age or max time horizon
                while(simulation.Agent.alive && simulation.Agent.age <= max_age && simulation.Agent.cal_year <= max_time_horizon)

                    # I have other event processes

                    # death - the last process
                    if simulation.Death.process(simulation.Agent)
                        # dead
                        simulation.Agent.alive = false
                        event_matrix[simulation.Agent.cal_year,simulation.Agent.age+1,simulation.Agent.sex+1,1] += 1
                    else
                        # still alive!
                        event_matrix[simulation.Agent.cal_year,simulation.Agent.age+1,simulation.Agent.sex+1,2] += 1
                        simulation.Agent.age += 1
                        simulation.Agent.cal_year += 1
                    end

                end # end while loop

            end # end for loop: agents

        end # end of begin timeit

    end # end for loop: cal year
    print_timer(to::TimerOutput)

    # convert the event matrix into a dictionary
    tmp_event_dict = Dict()
    for jj in 1:length(event_list)
        tmp_event_dict[event_list[jj]] = [event_matrix[:,:,1,jj],event_matrix[:,:,2,jj]]
    end

    # reshape event matrix into a dictionry of list of matrices
    simulation.OutcomeMatrix = Dict("n" => n_list,
                                    "outcome_matrix" => tmp_event_dict)

    print("\n Simulation finished. Check your simulation object for results.")
    return nothing
end

# test run
simulation = Simulation(111,2015,25,5000, agent, birth, death,Dict());

run(simulation)

You can directly download the files from: https://www.dropbox.com/sh/nrip5dc2rus6oto/AADIlZsngwjuUOwir9O2AeMsa?dl=0

I have just provided two modules, and my run function is very simple but you can imagine I would have many more processes and hence a complicated run function.

Besides any feedback you may have, I would like to know whether there is a better (in terms of readability, scalability, and efficiency) frameworks/tools I could use (e.g., software-wise, structure of "modules", julia types, built-in functions, etc).

If anything is unclear please let me know. Thanks!

Answer 1

That was a fun one. I only comment on the Julia parts, not the content, since the former was intense enough.

Deconstruction

abstract type Simulation_Module end

In Julia, PascalCase, as in SimulationModule, is preferred for types. For variable/function names, opinions differ between flatcase and snake_case.

---

mutable struct Agent{Z<:Bool,Y<:Int64}  <: Agent_Module
    sex::Z
    age::Y
    birth_year::Y
    cal_year::Y
    alive::Z
end

Those type parameters make no sense at all. Bool and Int64 are concrete types, so just use them directly: sex::Bool etc. I'd spell out calendar_year explicitely. Replace Int64 by Int unless you have specific reasons for using a certain size (Int is a uses the platforms word size).

Also, do think twice whether your types need to be mutable. Perhaps the algorithm works just as well with immutable types and the respective style.

---

function set_agent!(ag, sex, age, birth_year,cal_year,alive)
    ag.sex = sex;
    ag.age= age;
    ag.birth_year = birth_year;
    ag.cal_year = cal_year;
    ag.alive = alive;
    nothing
end

That's a weird one, but I understand it follows from making Agent mutable. I suggest to replace usages of set_agent! by updates with a new agent: simulation.agent = Agent(...).

If you insist on this function: take care of proper spacing (ag.age = age), and return the updated instance from it (return ag). Throw away those colons ;, too.

BTW, you load Setfield below. This function would be the ideal place to refactor using it.

---

mutable struct Birth <: Birth_Module
    parameters::Union{AbstractDict,Nothing}
    process::Function
end

Here we have the opposite case from Agent: AbstractDict and Function are abstract types, so you do want to make them parameters. Besides, abstracting the dict type might be an overkill; isn't there more specific information you have?

Furthermore, is it really needed to distinguish between nothing and an empty dictionary? Also, my point about mutability from above.

Suggestion:

mutable struct Birth{F, P} <: BirthModule
    parameters::Dict{Symbol, P}
    process::F
end

---

parameter_names = nothing

Uh. You know what you are doing there? You declare a mutable global variable, which is bad for type stability to begin with. And then you write to it multiple times from the top level, in every include!

Make the variable local wherever you use it.

---

birth = Birth(Dict_initializer(parameter_names),birth_process)

Same here. Make it local.

---

function Dict_initializer(parameter_names::Union{Nothing,Vector{Symbol}})
    isnothing(parameter_names) ? nothing : Dict(parameter_names .=> missing)
end

Style-wise, call it init_parameters or something.

And I suggest to get rid of it completely. It's ugly and will lead to weird types. Probably you want to avoid checking for the parameters being present in the simulation functions; you could refactor that by using convenience functions such as get! wherever the parameters are modified.

Or even better, implement immutable parameter updates using NamedTuple and SetField, with an appropriate initialization.

---

    Agent::Agent_Module
    Birth::Birth_Module
    Death::Death_Module

Fields should be in lower case.

In the case of Simulation, a mutable type seem reasonable to me. But perhaps a more elegant design is possible by separating simulation parameters (agent, time horizon, max age) that stays constant from a mutable "simulation state" type.

---

OutcomeMatrix::AbstractDict # needs to be a dictionary

And why is it called a matrix, then? Also, either make the dictionary type a parameter ({D<:AbstractDict}), or switch to a concrete type.

---

max_age::Int64 = simulation.max_age

The type assertion is not necessary. Also, you have SetField to do this for you!

---

simulation.n = (isnothing(simulation.n) ? 100 : simulation.n)

simulation.n is an Int64, and thus cannot be nothing.

Pass the initial n as a parameter to the run function, or put it into the simulation parameter type I proposed before, or at least make the number a global const.

---

max_time_horizon::Int64 = simulation.time_horizon

Use consistent naming here.

---

cal_years::Vector{Int64} = collect(min_cal_year:1:max_cal_year)

collect(min_cal_year:max_cal_year) should suffice. But think about whether you need to allocate this vector in the first place. Aren't you only iterating over it? Then the range is enough. Or could you even just recalculate it in every step from the initial year and the current n?

---

n_list = zeros(Int64,simulation.time_horizon,2)

Not a good name. It's neither a list (but a Vector), nor has it an appearent connection to n.

Also, it has fixed size. You could replace it by a tuple or two variables.

Watch your spaces!

---

event_list::Vector{String} = ["death","alive"]

Symbols would be preferred: [:death, :alive].

---

    tmp_n::Int64 = 0
    tmp_cal_year_index::Int64 = 0

Better names for that, please! Or get rid of the variable all together.

---

event_matrix = zeros(Int64,length(cal_years),max_age,2,length(event_list))

Maybe this can get a better name, too... events? And it seems that you didn't completely decide whether to make the kinds of events fixed or parametrizable. If there are only ever going to be two kinds, use a const N_EVENTS = 2. Otherwise, parametrize the simulation by event types (and then change how you handle event_list).

Also, that 2 needs to be a constant too, I have no idea what it the dimension stands for. Perhaps use a type like AxisArray to make all dimensions more explicit.

Spacing, again!

---

n_list[tmp_cal_year_index,simulation.Agent.sex+1] +=1

Ah, so this should have been called agent_counts or like that. Again, spacing!

n_list[tmp_cal_year_index, simulation.agent.sex + 1] += 1

And wait, what? The sex is a Bool. Why are you incrementing it? Make the semantics of using a Bool to select the column more explicit. Perhaps there is an array type that can do that. You could create your own even (or give up the social construction of there being a limited number of two distinct sexes :))

---

if simulation.Death.process(simulation.Agent)

This looks like an updating operation. In that case, the name of the field of the death module should be process!.

This also shows that you should add more docstrings. Document the interface that process! needs to fulfill somewhere, too!

---

# convert the event matrix into a dictionary
    tmp_event_dict = Dict()
    for jj in 1:length(event_list)
        tmp_event_dict[event_list[jj]] = [event_matrix[:,:,1,jj],event_matrix[:,:,2,jj]]
    end

    # reshape event matrix into a dictionry of list of matrices
    simulation.OutcomeMatrix = Dict("n" => n_list,
                                    "outcome_matrix" => tmp_event_dict)

This looks a bit gross. Suggestion: have a separate type for simulation results. Then drop OutcomeMatrix and instead of writing the simulation result to the simulation object, return a simulation result object.

Apart from that:

• tmp_ is a really bad prefix. Try to construct the thing directly:

  outcome_matrix = [(event_matrix[:,:,1,jj], event_matrix[:,:,2,jj]) for jj in eachindex(event_list)]
  

  This also uses tuples instead of constant-length arrays, and avoids an untyped Dict().

  Additionally, you copy over the whole event matrix into slices. That's bad for memory. You could use views, or make the logic part of the suggested result type -- pack the event matrix into it, and add some nice custom accessors instead of the dictionary conversion.

• Once more, OutcomeMatrix is obviously not a matrix. If you do something like this at all, use Symbols in the dict. Even better, use a NamedTuple instead: (;n = n_list, outcome_matrix = outcome_matrix), or in newer Julias (;n = n_list, outcome_matrix) with the naming inferred automatically.

---

simulation = Simulation(111,2015,25,5000, agent, birth, death,Dict());

Again, this looks like Fortran with more indirection: you pass in a mutable object to use as the output. Make the simulation result the return value of run.

Remove the colon.

More general design suggestions

1. Think about whether you need so much abstraction. E.g., the abstract base types for everything.
2. Perhaps try to rewrite run as a combination of an init! and a step! function on a mutable SimulationState object. Then, derive an iterator from it so that you can do for step in run(simulation) ....