Julia ❤️ ABM: Issue 3

What a difference a day makes

In the previous issue we have laid out the basic events that happen in our model - both as verbal descriptions in plain English and as algorithmic representations in simple Julia code. Sally and Richard are now not merely static pieces of data anymore but they can act in three different ways: work, eat, and trade.

However, today's lesson content will be a bit more technical in nature and deal with how we work with Julia. So far we had just written out the pieces of our modelling code in long blocks of for-loops that stood on their own and weren't really connected to each other. While it works great for simple examples, this will become increasingly harder to read. Furthermore, it is relatively complicated to reason about which pieces of our code belong to which category of actions. To take care of this, today we will primarily learn about two things:

1. Writing our code in plain text files and including them in a running Julia REPL.
2. Encapsulating our code in functions takes care of this and provides us with more options to easily structure our model code and even improve its performance.

So let us not hesitate any longer and start the refactoring.

Storing code in text files

In the previous issues, we simply had a Julia REPL running and directly input all the interesting bits of code. And every time we closed the REPL at the end of a day or whenever we feel like we had learned enough about agent-based modelling in Julia, we had to re-enter all the previous code blocks from earlier sessions. As you can imagine, this quickly becomes a time sink of unexpected dimensions but luckily there is an easy way around this.

Whenever we have a piece of code that we want to reuse every time we work on our ABM, we should store said code in plain text in its own file. It's always a good idea to organise your projects into separate folders, so we will now start to do that by using built-in functions inside the Julia REPL.

First we start our Julia REPL again by calling julia in the command line. We then create a new directory called julia-loves-abm (or any other name for that matter) by using the mkdir function (read: make directory).

julia> mkdir("julia-loves-abm")  
"julia-loves-abm"

Next we switch our current working directory to that newly created folder with the cd function (read: change directory)...

julia> cd("julia-loves-abm/")

... and make sure everything worked as expected by checking that our current working directory was properly updated with the pwd function (read: print working directory). Note that the output of the pwd() function will likely look a bit different for you and will depend on the folder from which you started Julia in as well as your operating system (Linux, Mac, Windows, etc.). Personally, I like to put all my coding projects into a Code/ folder so that I can quickly find all my old and new projects. Anyways, the important part about it is that the last bit of the returned string says julia-loves-abm which is the newly created folder.

julia> pwd()  
"/home/fbanning/Code/julia-loves-abm"

Now that we have created and switched to our dedicated project folder, we can create a new text file with the touch function.

julia> touch("WorkEatTradeRepeat.jl")  
"WorkEatTradeRepeat.jl"

While the name for the file is mostly arbitrary, you might have noticed that we've also added a specific file extension, namely .jl. This is the common identifier for text files containing Julia code.

Alright, we have successfully created a new file that we can open with the text editor of our choice (e.g. notepad++, Kate, VS Code, or any other that you prefer). Julia also provides a built-in way to open your preferred text editor from inside the REPL with the edit function which even provides some nice extra functionality to open a file on a given line or at the beginning of a function definition. However, I've found that it does not always work straight out of the box and might require some manual setup first. Feel free to try it out if you're curious. :)

So far the file is completely empty and devoid of any meaningful content. We can go ahead and copy/paste most of our code from the previous issue. Here's a quick recap how that might look:

mutable struct Agent
    job::String
    pantry_fish::Float64
    pantry_bread::Float64
    hunger_fish::Float64
    hunger_bread::Float64
end

sally = Agent( "Baker", 10, 10, 0, 0)
richard = Agent("Fisher", 10, 10, 0, 0)
agents = [sally, richard]

for agent in agents
    if agent.job == "Baker"
        fresh_breads = rand(1:20)
        agent.pantry_bread += fresh_breads
        baking_hours = LinRange(8, 10, 20)[fresh_breads]
    elseif agent.job == "Fisher"
        fishing_hours = rand(1:5)
        fresh_fish = rand([0, 0, 0, 0, 1, 1, 1, 2, 2, 3], fishing_hours)
        agent.pantry_fish += sum(fresh_fish)
        other_work_hours = rand(1:3)
    end
    # get hungry from working
    agent.hunger_bread += 1.0
    agent.hunger_fish += 0.5
end

for agent in agents
    # completely satisfy hunger when there's enough food in pantry
    if agent.pantry_fish >= agent.hunger_fish
        agent.pantry_fish -= agent.hunger_fish
        agent.hunger_fish = 0
    end
    if agent.pantry_bread >= agent.hunger_bread
        agent.pantry_bread -= agent.hunger_bread
        agent.hunger_bread = 0
    end
    
    # partially satisfy hunger until pantry is empty
    if 0 < agent.pantry_fish <= agent.hunger_fish
        agent.hunger_fish -= agent.pantry_fish
        agent.pantry_fish = 0
    end
    if 0 < agent.pantry_bread <= agent.hunger_bread
        agent.hunger_bread -= agent.pantry_bread
        agent.pantry_bread = 0
    end
end

if sally.pantry_fish < 1.5
    sally_wants_fish = 1.5 - sally.pantry_fish
    traded_fish = min(sally_wants_fish, sally.pantry_bread / 2, richard.pantry_fish)
    traded_bread = 2 * traded_fish

    sally.pantry_bread -= traded_bread
    sally.pantry_fish += traded_fish
    richard.pantry_fish -= traded_fish
    richard.pantry_bread += traded_bread
end

if richard.pantry_bread < 3.0
    richard_wants_bread = 3.0 - richard.pantry_bread
    traded_bread = min(richard_wants_bread, sally.pantry_bread, richard.pantry_fish * 2)
    traded_fish = traded_bread / 2

    sally.pantry_bread -= traded_bread
    sally.pantry_fish += traded_fish
    richard.pantry_fish -= traded_fish
    richard.pantry_bread += traded_bread
end

Go ahead and save this file now.

Back in the Julia REPL we can call the include function to load our .jl-file and use its content in the running process.

julia> include("WorkEatTradeRepeat.jl")

Everything from inside the file has now been read and evaluated by Julia from top to bottom. What happened here is fundamentally the same as what happened before when we directly typed things into the REPL:

1. Define the agent variables,
2. create Sally and Richard and store them in a vector, and finally
3. let them work, eat, and trade once.

Note that this also means that both Sally and Richard continue to exist after we have included our file. If we now check on their variables, we can see their current states as we've been used to from before.

julia> richard
Agent("Fisher", 19.0, 10.0, 0.0, 0.0)  

julia> sally
Agent("Baker", 10.0, 30.0, 0.0, 0.0)

Now whenever we change something in our WorkEatTradeRepeat.jl file, we have to repeat the steps above: First save it and then include it in the Julia REPL for its content to be evaluated again.

Encapsulating code in functions

Currently we have our logic wrapped in for agent in agents loops that sequentially executes the same code for each of our agents. We now separate the model logic from the repetitions by refactoring parts of it into their own functions.

Let's start with putting the code that describes Sally's and Richard's work days into its own function called work!. It takes any agent as its sole function argument, so it doesn't matter if we pass sally or richard to it.

function work!(agent)
    if agent.job == "Baker"
        fresh_breads = rand(1:20)
        agent.pantry_bread += fresh_breads
        baking_hours = LinRange(8, 10, 20)[fresh_breads]
        work_hours = baking_hours
    elseif agent.job == "Fisher"
        fishing_hours = rand(1:5)
        fresh_fish = rand([0, 0, 0, 0, 1, 1, 1, 2, 2, 3], fishing_hours)
        agent.pantry_fish += sum(fresh_fish)
        other_work_hours = rand(1:3)
        work_hours = fishing_hours + other_work_hours
    end
    # get hungry from working
    agent.hunger_bread += 1.0 + work_hours / 8
    agent.hunger_fish += 0.5 + work_hours / 16

    return agent
end

There's also a bit of new functionality in there from now on: The more hours Sally and Richard spend working, the more calories they will burn and therefore also get hungrier. There is a basic need for bread and fish every day, even if they would not work at all, and their hunger for bread and fish increases with different intensity according to their daily workload.

After making this function available in the Julia REPL (i.e. save the file and include it), it is now available for us to use it repeatedly. To recreate our code from the last issue, we can now simply write

for agent in agents
    work!(agent)
end

to achieve exactly the same outcome as before. This makes it a lot clearer what is happening inside the for-loop at a quick glance. And while it may seem like a very minor change at first, breaking up the code of our model into smaller pieces is a fast and easy way to maintain its readability. In this case, it allows us to jump directly to the function that defines the working behaviour of our agents and modify just that small part we're currently concerned with.

Needless to say, we will do the same with the code blocks that let our agents eat! and trade!.

function eat!(agent)
    # completely satisfy hunger when there's enough food in pantry
    if agent.pantry_fish >= agent.hunger_fish
        agent.pantry_fish -= agent.hunger_fish
        agent.hunger_fish = 0
    end
    if agent.pantry_bread >= agent.hunger_bread
        agent.pantry_bread -= agent.hunger_bread
        agent.hunger_bread = 0
    end

    # partially satisfy hunger until pantry is empty
    if 0 < agent.pantry_fish <= agent.hunger_fish
        agent.hunger_fish -= agent.pantry_fish
        agent.pantry_fish = 0
    end
    if 0 < agent.pantry_bread <= agent.hunger_bread
        agent.hunger_bread -= agent.pantry_bread
        agent.pantry_bread = 0
    end

    return agent
end

For the trade! function we do not use the sally and richard variables anymore like we previously did. Instead we opt for a more general description of the trading agents by calling the function arguments agent1 and agent2. Copying and pasting what we have coded in the previous issue, we can simply replace any occurrence of sally with agent1 and every richard with agent2. You can think of this as creating aliases for the two agents in question which puts our code in a more generalised form and allows us to later reuse it for agents that are different from our two friends.⊕ Ooohh, such foreshadowing, much suspense!

function trade!(agent1, agent2)
    if agent1.pantry_fish < 1.5
        agent1_wants_fish = 1.5 - agent1.pantry_fish
        traded_fish = min(agent1_wants_fish, agent1.pantry_bread / 2, agent2.pantry_fish)
        traded_bread = 2 * traded_fish

        agent1.pantry_bread -= traded_bread
        agent1.pantry_fish += traded_fish
        agent2.pantry_fish -= traded_fish
        agent2.pantry_bread += traded_bread
    end

    if agent2.pantry_bread < 3.0
        agent2_wants_bread = 3.0 - agent2.pantry_bread
        traded_bread = min(agent2_wants_bread, agent1.pantry_bread, agent2.pantry_fish * 2)
        traded_fish = traded_bread / 2

        agent1.pantry_bread -= traded_bread
        agent1.pantry_fish += traded_fish
        agent2.pantry_fish -= traded_fish
        agent2.pantry_bread += traded_bread
    end

    return agent1, agent2
end

The following subsections explain some details about how Julia functions work and how they are most commonly shaped. This is not at all exhaustive with respect to the things that can be done but instead highlights a few things that agent-based modellers might encounter more often than not. If you are interested to delve a bit deeper into functions, you can have a look at the section in the official Julia docs as usual. Should you feel like skipping these details or just already know most of them, go ahead and continue reading here.

In-place functions

You may have noticed that all the function names that we have written so far end with an exclamation mark ! (also sometimes referred to as a "bang"). This was a deliberate decision as it follows a common convention in the Julia community to mark in-place functions. These functions mutate the values of their arguments instead of first making a copy of them or creating new variables altogether.

Using the same Agent struct form above, let's illustrate this in a simple code example with a very hungry agent:

julia> a = Agent("Scientist", 0, 0, 9000, 9000)  
Agent("Scientist", 0.0, 0.0, 9000.0, 9000.0)

julia> function add_to_hunger(a, b)
           a.hunger_fish + b
           a.hunger_bread + b
           return a
       end
add_to_hunger (generic function with 1 method)

julia> add_to_hunger(a, 1)

julia> a
Agent("Scientist", 0.0, 0.0, 9000.0, 9000.0)

julia> function add_to_hunger!(a, b)
           a.hunger_fish += b
           a.hunger_bread += b
           return a
       end
add_to_hunger! (generic function with 1 method)

julia> add_to_hunger!(a, 1)
Agent("Scientist", 0.0, 0.0, 9001.0, 9001.0)

julia> a
Agent("Scientist", 0.0, 0.0, 9001.0, 9001.0)

As you can see, the exclamation mark in the function name has actually no effect on what the function itself does to its arguments. They are even named almost the same, except for the ! at the end of the function name. However, add_to_hunger did not change the outside variable a while add_to_hunger! did change the value of a.

• The add_to_hunger function only took the values of the two hunger fields of our starving scientist and added b to them. But in fact the agent's hunger variables were not updated in the process and the newly calculated numbers were dropped after finishing the function execution. At the end, the function returned the a variable in the same state as it was before.
• The add_to_hunger! function assigned the resulting values (9001.0 in both cases) to the agent fields and by that also updated the variable a that was initially defined outside of the function body. And although not strictly necessary, it is indeed considered best-practice to have functions that modify its arguments in-place return the modified variables.

Even though the function names may have hinted at this behaviour, we cannot always rely on the function behaving as expected when we see a ! at the end of a function name. It is solely a convention widely used throughout the Julia community. Also note that the arguments have to be of a mutable type, just like Agent in our example above.

In general it can be said that we will most often deal with in-place functions in agent-based modelling since we already have well-defined existing agent variables and model properties that will be modified over time (e.g. at each time step in our simulation). Still, it's good to know about the distinction between them, especially when making use of code that was written by others.

Specialised function methods

When writing a function, we can in fact define multiple specialised versions of it which share the same name. These are called methods and their creation can happen in two distinct ways: Either by creating the function with the same (keyword) arguments and annotating them with a specific type. Or by adding more (keyword) arguments altogether.

Option 1 can be demonstrated in the ABM context when we think about multiple types of agents. So far we have distinguished between Sally the baker and Richard the fisherman by checking for their current value in the job field of the Agent struct. Inside the work! function we have a simple conditional check on agent.job and then run the appropriate commands for that work. This works just fine in our case and doesn't require a lot of computational resources to find out which job the current agent has.

However, we could have also created specialised agent types to reflect different jobs and then let Julia dispatch the agent to the correct function method. This may sound a bit confusing at first, so let's write it out in code to explain the concept behind it.

First we define two specialised structs for Bakers and Fishers.

julia> mutable struct Baker  
          pantry_fish::Float64  
          pantry_bread::Float64  
          hunger_fish::Float64  
          hunger_bread::Float64  
      end  
  
julia> mutable struct Fisher  
          pantry_fish::Float64  
          pantry_bread::Float64  
          hunger_fish::Float64  
          hunger_bread::Float64  
      end

As you can see, the struct fields are exactly the same but their names differ. Now when we create Sally the baker and Richard the Fisher again, we can use the custom types Baker and Fisher and have essentially the same information as before when we were using the job field in the Agent type.

julia> sally = Baker(10, 10, 0, 0)
Baker(10.0, 10.0, 0.0, 0.0)

julia> richard = Fisher(10, 10, 0, 0)
Baker(10.0, 10.0, 0.0, 0.0)

One neat side-effect of this is, that we can now use some built-in syntax to check which job each agent has.

julia> typeof(sally)
Baker

julia> sally isa Baker
true

Next we copy the work! function from above but split it into two different methods - one for the work of a Baker and one for the work of a Fisher.

julia> function work!(agent::Baker)
           fresh_breads = rand(1:20)
           agent.pantry_bread += fresh_breads
           baking_hours = LinRange(8, 10, 20)[fresh_breads]
       
           return agent
       end
work! (generic function with 2 methods)
  
julia> function work!(agent::Fisher)
           fishing_hours = rand(1:5)
           fresh_fish = rand([0, 0, 0, 0, 1, 1, 1, 2, 2, 3], fishing_hours)
           agent.pantry_fish += sum(fresh_fish)
           other_work_hours = rand(1:3)
       
           return agent
       end
work! (generic function with 3 methods)

Note how we have essentially written the same functions as before but annotated a type for the single argument that the function takes. This is done via the ::T syntax behind the argument where T can be any custom (Baker) or built-in (Vector) type. The last line tells us that work! now is defined as a function with 3 methods, hinting at the fact that we have indeed succeeded in what we set out to achieve.

Now comes the amazing part of this setup: We can call the work! function by passing an instance of our specialised agent types to it and Julia will automatically pick the correct method.

julia> work!(sally)  
Baker(10.0, 19.0, 0.0, 0.0)

julia> work!(richard)  
Fisher(11.0, 10.0, 0.0, 0.0)

Hey, look at that! Sally has baked nine loaves of bread and Richard has caught one fish that day. This makes things very easy for us and we can now create specialised methods for agent behaviour functions throughout our model as long as the behaviour is based on the agents' jobs.

A day in the life

Wow, that was quite a lot of information, wasn't it? We have restructured our agent types by the jobs (Baker and Fisher) and now understand the difference between mutating (denoted with ! at the end of the function name) and non-mutating functions. But we're not quite finished yet with today's issue, so let's push on to bring it all together.

We can bundle the previously defined work!, eat!, and trade! functions into one function that will advance our simulation by one time step. In our case this equals to one day. Working and eating happen independently from other agents which allows us to loop over all our agents one after another.

Trading, however, happens between two agents which is also why the trade! function takes two agents as its arguments. So instead of iterating over our two agents one by one, we call the trade! function once and give it a pair of elements from the agents vector to work with. Currently this vector only consists of Sally and Richard, therefore we can retrieve these agents by directly indexing into the vector with the square bracket notation []. For example, agents[1] gives us the first element of that vector.

function a_day_in_the_life!(agents)
    for agent in agents
        work!(agent)
        eat!(agent)
    end
    
    trade!(agents[1], agents[2])
    
    return agents
end

To simulate 30 days in a row, we can simply call this function repeatedly after creating our two agents Sally and Richard.

richard = Fisher(10, 10, 0, 0)
sally   = Baker(10, 10, 0, 0)
agents  = [sally, richard]
for day in 1:30
    a_day_in_the_life!(agents)
end

Hey, the building blocks are slowly coming together and the code in our WorkTradeEatRepeat.jl file is already starting to look like an actual agent-based model.

Note that the outer for-loop iterates over the days/steps/ticks of our simulation (for day in 1:30). The for-loop inside the a_day_in_the_life! function then does exactly the working and eating processes that we've described above for each agent in agents in a sequential manner (i.e. those are independent processes where the results do not rely on the state of other agent variables or the model environment). Then, at last during that day, we execute the trading code block. Once all of this has happened, our program goes back to the beginning of the outer for-loop and continuously advances our simulation state for the remaining 29 days.

Trading every third day

But wait, while our agents work and eat on every day, it's only on every third day that they also have the opportunity to trade. To find out which of the potentially hundreds of days that we want to simulate is a "third day", we can use the modulo operator % in Julia. It provides us with the remainder of a mathematical division. Here are a few examples of its usage with comments to describe the results:

julia> 3 % 1 # special case, always returns 0 because everything is a multiple of 1
0

julia> 3 % 2 # 3 modulo 2 has a quotient of 1 and a remainder of 1 which is returned
1

julia> 3 % 3 # returns 0 because dividend is a multiple of divisor
0

julia> 3 % 4 # returns dividend when divisor is bigger
3

If the current day is a trading day (i.e. a day that is a multiple of 3), then that remainder will be equal to zero. Thus, we can use this result for a logical comparison in our code and only execute the block of code related to trading on the correct days. Incorporating this knowledge in our a_day_in_the_life! function looks like this:

function a_day_in_the_life!(agents, day)
    for agent in agents
        work!(agent)
        eat!(agent)
    end

    if day % 3 == 0
        trade!(agents[1], agents[2])
    end
end

Note that we had to include a second keyword day in the function definition. This is necessary because we can generally only access those variables inside the function that have been passed to it when it was called. While it may seem annoying at first, this is indeed a great feature and not (only) a nuisance. If you're interested in the issue of scoping and namespaces, you're highly encouraged to wager a deep dive into the Julia documentation again. As always, it might be helpful to know these things but for the purpose of being able to follow along with this tutorial series, it is not particularly necessary to understand more than the basics behind it. And knowing about the current day is indeed crucial for our model logic to work as intended.

In the above-mentioned case of a 30 day simulation, our main loop has also changed a little bit as a consequence because we need to pass the current day as the second argument to the function call.

richard = Fisher(10, 10, 0, 0)
sally   = Baker(10, 10, 0, 0)
agents  = [sally, richard]
for day in 1:30
    a_day_in_the_life!(agents, day)
end

Summary and outlook

In this issue of the "Julia ❤️ ABM" series, we have learned about...

• how to make our model code more readable and increase its maintainability by encapsulating code blocks in separate functions,
• leaving the primarily REPL-based workflow behind and adding our code to locally saved .jl-files,
• the modulo operator %,
• bundling the model logic in a stepping function.

All of this is pretty great because it allows us to separate the model's sub-processes from the broader picture of how our model advances over time. This is very similar to the separation of intermediary steps inside NetLogo's go procedure into their own custom procedures. Reading the code we wrote over the course of this issue, we can now see at a glance:

1. What is happening in one simulated day of our model (thanks to carefully chosen telling function names), and
2. in which order things are happening in our model (thanks to the separation with conditional logic).

For convenience, here's how the WorkEatTradeRepeat.jl file might look like up until this point.

In the next issue, we will deviate a bit from the initial concept of laying out the functional framework of our agent-based model. Instead our focus will continue to shift towards improving our coding environment. We will think about how to generalise the initialisation of our model, expand its scope from just two agents to a whole lot more, make use of built-in and publicly available packages to make things as easy as possible for us, and finally how to collect data from our simulation runs.