Pure Functional Epidemics An Agent-Based Approach

Agent-Based Simulation (ABS) is a methodology in which a system is simulated in a bottom-up approach by modelling the micro interactions of its constituting parts, called agents, out of which the global system behaviour emerges. So far mainly object-oriented techniques and languages have been used in ABS. Using the SIR model of epidemiology, which simulates the spreading of an infectious disease through a population, we demonstrate how to use pure Functional Reactive Programming to implement ABS. With our approach we can guarantee the reproducibility of the simulation at compile time and rule out specific classes of run-time bugs, something that is not possible with traditional object-oriented languages. Also, we found that the representation in a purely functional format is conceptually quite elegant and opens the way to formally reason about ABS.


INTRODUCTION
The traditional approach to Agent-Based Simulation (ABS) has so far always been object-oriented techniques, due to the influence of the seminal work of Epstein et al [15] in which the authors claim "[..] object-oriented programming to be a particularly natural development environment for Sugarscape specifically and artificial societies generally [..]" (p.179).This work established the metaphor in the ABS community, that agents map naturally to objects [31], which still holds up today.
In this paper we challenge this metaphor and explore ways of approaching ABS in a pure (lack of implicit side-effects) functional way using Haskell.By doing this we expect to leverage the benefits of pure functional programming [21]: higher expressivity through declarative code, being polymorph and explicit about side-effects through monads, more robust and less susceptible to bugs due to explicit data flow and lack of implicit side-effects.
As use case we introduce the SIR model of epidemiology with which one can simulate epidemics, that is the spreading of an infectious disease through a population, in a realistic way.
Over the course of three steps, we derive all necessary concepts required for a full agent-based implementation.We start with a Functional Reactive Programming (FRP) [46] solution using Yampa [20] to introduce most of the general concepts and then make the transition to Monadic Stream Functions (MSF) [34] which allow us to add more advanced concepts of ABS to pure functional programming.
The aim of this paper is to show how ABS can be implemented in pure Haskell and what the benefits and drawbacks are.By doing this we give the reader a good understanding of what ABS is, what the challenges are when implementing it and how we solve these in our approach.
The contributions of this paper are: • We present an approach to ABS using declarative analysis with FRP in which we systematically introduce the concepts of ABS to pure functional programming in a step-by-step approach.Also this work presents a new field of application to FRP as to the best of our knowledge the application of FRP to ABS (on a technical level) has not been addressed before.The result of using FRP allows expressing continuous time-semantics in a very clear, compositional and declarative way, abstracting away the low-level details of time-stepping and progress of time within an agent.• Our approach can guarantee reproducibility already at compile time, which means that repeated runs of the simulation with the same initial conditions will always result in the same dynamics, something highly desirable in simulation in general.This can only be achieved through purity, which guarantees the absence of implicit side-effects, which allows to rule out non-deterministic influences at compile time through the strong static type system, something not possible with traditional object-oriented approaches.Further, through purity and the strong static type system, we can rule out important classes of run-time bugs e.g.related to dynamic typing, and the lack of implicit data-dependencies which are common in traditional imperative object-oriented approaches.
• Using pure functional programming, we can enforce the correct semantics of agent execution through types where we demonstrate that this allows us to have both, sequential monadic behaviour, and agents acting conceptually at the same time in lock-step, something not possible using traditional object-oriented approaches.
In Section 2 we define Agent-Based Simulation, introduce Functional Reactive Programming, Arrowized programming and Monadic Stream Functions, because our approach builds heavily on these concepts.In Section 3 we introduce the SIR model of epidemiology as an example model to explain the concepts of ABS.The heart of the paper is Section 4 in which we derive the concepts of a pure functional approach to ABS in three steps, using the SIR model.Section 5 discusses related work.Finally, we draw conclusions and discuss issues in Section 6 and point to further research in Section 7.

BACKGROUND 2.1 Agent-Based Simulation
Agent-Based Simulation is a methodology to model and simulate a system where the global behaviour may be unknown but the behaviour and interactions of the parts making up the system is known.Those parts, called agents, are modelled and simulated, out of which then the aggregate global behaviour of the whole system emerges.
So, the central aspect of ABS is the concept of an agent which can be understood as a metaphor for a pro-active unit, situated in an environment, able to spawn new agents and interacting with other agents in some neighbourhood by exchange of messages.
We informally assume the following about our agents [26,39,47]: • They are uniquely addressable entities with some internal state over which they have full, exclusive control.• They are pro-active, which means they can initiate actions on their own e.g.change their internal state, send messages, create new agents, terminate themselves.• They are situated in an environment and can interact with it.
• They can interact with other agents situated in the same environment by means of messaging.
Epstein [14] identifies ABS to be especially applicable for analysing "spatially distributed systems of heterogeneous autonomous actors with bounded information and computing capacity".They exhibit the following properties: • Linearity & Non-Linearity -actions of agents can lead to non-linear behaviour of the system.• Time -agents act over time, which is also the source of their pro-activity.

Functional Reactive Programming
Functional Reactive Programming is a way to implement systems with continuous and discrete time-semantics in pure functional languages.There are many different approaches and implementations but in our approach we use Arrowized FRP [22,23] as implemented in the library Yampa [9,20,29].The central concept in Arrowized FRP is the Signal Function (SF), which can be understood as a process over time which maps an input-to an output-signal.A signal can be understood as a value which varies over time.Thus, signal functions have an awareness of the passing of time by having access to ∆t which are positive time-steps, the system is sampled with.
Yampa provides a number of combinators for expressing timesemantics, events and state-changes of the system.They allow to change system behaviour in case of events, run signal functions and generate stochastic events and random-number streams.We shortly discuss the relevant combinators and concepts we use throughout the paper.For a more in-depth discussion we refer to [9,20,29].
Event.An event in FRP is an occurrence at a specific point in time, which has no duration e.g. the recovery of an infected agent.Yampa represents events through the Event type, which is programmatically equivalent to the Maybe type.
Dynamic behaviour.To change the behaviour of a signal function at an occurrence of an event during run-time, (amongst others) the combinator switch :: SF a (b, Event c) → (c → SF a b) → SF a b is provided.It takes a signal function, which is run until it generates an event.When this event occurs, the function in the second argument is evaluated, which receives the data of the event and has to return the new signal function, which will then replace the previous one.Note that the semantics of switch are that the signal function, into which is switched, is also executed at the time of switching.
Randomness.In ABS, often there is the need to generate stochastic events, which occur based on e.g. an exponential distribution.
Yampa provides the combinator occasionally :: RandomGen g ⇒ g → Time → b → SF a (Event b) for this.It takes a random-number generator, a rate and a value the stochastic event will carry.It generates events on average with the given rate.Note that at most one event will be generated and no 'backlog' is kept.This means that when this function is not sampled with a sufficiently high frequency, depending on the rate, it will lose events.
Yampa also provides the combinator noise :: (RandomGen g, Random b) ⇒ g → SF a b, which generates a stream of noise by returning a random number in the default range for the type b.
Running signal functions.To purely run a signal function, Yampa provides the function embed :: SF a b → (a, [(DTime, Maybe a)]) → [b], which allows to run an SF for a given number of steps where in each step one provides the ∆t and an input a.The function then returns the output of the signal function for each step.Note that the input is optional, indicated by Maybe.In the first step at t = 0, the initial a is applied and whenever the input is Nothing in subsequent steps, the most recent a is re-used.

Arrowized programming
Yampa's signal functions are arrows, requiring us to program with arrows.Arrows are a generalisation of monads, which in addition to the already familiar parameterisation over the output type, allow parameterisation over their input type as well [22,23].
In general, arrows can be understood to be computations that represent processes, which have an input of a specific type, process it and output a new type.This is the reason why Yampa is using arrows to represent their signal functions: the concept of processes, which signal functions are, maps naturally to arrows.
There exists a number of arrow combinators, which allow arrowized programing in a point-free style but due to lack of space we will not discuss them here.Instead we make use of Paterson's do-notation for arrows [32], which makes code more readable as it allows us to program with points.
To show how arrowized programming works, we implement a simple signal function, which calculates the acceleration of a falling mass on its vertical axis as an example [35].To create an arrow, the proc keyword is used, which binds a variable after which the do of Patersons do-notation [32] follows.Using the signal function integral :: SF v v of Yampa, which integrates the input value over time using the rectangle rule, we calculate the current velocity and the position based on the initial position p0 and velocity v0.The <<< is one of the arrow combinators, which composes two arrow computations and arr simply lifts a pure function into an arrow.To pass an input to an arrow, -< is used and <to bind the result of an arrow computation to a variable.Finally to return a value from an arrow, returnA is used.

Monadic Stream Functions
Monadic Stream Functions (MSF) are a generalisation of Yampa's signal functions with additional combinators to control and stack side effects.An MSF is a polymorphic type and an evaluation function, which applies an MSF to an input and returns an output and a continuation, both in a monadic context [34]: To show how arrowized programming with MSFs works, we extend the falling mass example from above to incorporate monads.In this example we assume that in each step we want to accelerate our velocity v not by the gravity constant anymore but by a random number in the range of 0 to 9.81.Further we want to count the number of steps it takes us to hit the floor, that is when position p is less than or equal 0. Also when hitting the floor we want to print a debug message to the console with the velocity by which the mass has hit the floor and how many steps it took.
We define a corresponding monad stack with IO as the innermost Monad, followed by a RandT transformer for drawing randomnumbers and finally a StateT transformer to count the number of steps we compute.We can access the monadic functions using arrM in case we need to pass an argument and _arrM in case no argument to the monadic function is needed: --drawing random number for our gravity range r <-arrM_ (lift $ lift $ getRandomR (0, 9.81)) -< () v <-arr (+v0) <<< integral -< (-r) p <-arr (+p0) <<< integral -< v --count steps arrM_ (lift (modify (+1))) -< () if p > 0 then returnA -< p --we have hit the floor else do --get number of steps s <-arrM_ (lift get) -< () --write to console arrM (liftIO .putStrLn) -< "hit floor with v " ++ show v ++ " after " ++ show s ++ " steps" returnA -< p To run the fallingMassMSF function until it hits the floor we proceed as follows:  extra lift for accessing StateT and RandT.Thus unMSF returns a computation in the ReaderT Double Monad, which we need to peel away using runReaderT.This then results in a StateT Int computation, which we evaluate by using runStateT and the current number of steps as state.This then results in another monadic computation of RandT Monad, which we evaluate using runRandT.This finally returns an IO computation, which we simply evaluate to arrive at the final result.

THE SIR MODEL
To explain the concepts of ABS and of our pure functional approach to it, we introduce the SIR model as a motivating example and use-case for our implementation.It is a very well studied and understood compartment model from epidemiology [25], which allows to simulate the dynamics of an infectious disease like influenza, tuberculosis, chicken pox, rubella and measles spreading through a population [13].
In this model, people in a population of size N can be in either one of three states Susceptible, Infected or Recovered at a particular time, where it is assumed that initially there is at least one infected person in the population.People interact with each other on average with a given rate of β per time-unit and become infected with a given probability γ when interacting with an infected person.When infected, a person recovers on average after δ time-units and is then immune to further infections.An interaction between infected persons does not lead to re-infection, thus these interactions are ignored in this model.This definition gives rise to three compartments with the transitions seen in Figure 1.
This model was also formalized using System Dynamics (SD) [36].In SD one models a system through differential equations, allowing to conveniently express continuous systems, which change over time, solving them by numerically integrating over time, which gives then rise to the dynamics.We won't go into detail here and provide the dynamics of such a solution for reference purposes, shown in Figure 2.

An Agent-Based approach
The approach of mapping the SIR model to an ABS is to discretize the population and model each person in the population as an individual agent.The transitions between the states are happening due to discrete events caused both by interactions amongst the agents and time-outs.The major advantage of ABS is that it allows to incorporate spatiality as shown in Section 4.3 and simulate heterogenity of population e.g.different sex, age.This is not possible with other simulation methods e.g.SD or Discrete Event Simulation [48].
According to the model, every agent makes on average contact with β random other agents per time unit.In ABS we can only contact discrete agents thus we model this by generating a random event on average every 1 β time units.We need to sample from an exponential distribution because the rate is proportional to the size of the population [5].Note that an agent does not know the other agents' state when making contact with it, thus we need a mechanism in which agents reveal their state in which they are in at the moment of making contact.This mechanism is an implementation detail, which we will derive in our implementation steps.For now we only assume that agents can make contact with each other somehow.

DERIVING A PURE FUNCTIONAL APPROACH
In [42] two fundamental problems of implementing an ABS from a programming-language agnostic point of view is discussed.The first problem is how agents can be pro-active and the second how interactions and communication between agents can happen.For agents to be pro-active, they must be able to perceive the passing of time, which means there must be a concept of an agent-process, which executes over time.Interactions between agents can be reduced to the problem of how an agent can expose information about its internal state which can be perceived by other agents.
Both problems are strongly related to the semantics of a model and the authors show that it is of fundamental importance to match the update-strategy with the semantics of the model -the order in which agents are updated and actions of agents are visible can make a big difference and need to match the model semantics.The authors identify four different update-strategies, of which the parallel update-strategy matches the semantics of the agent-based SIR model due to the underlying roots in the System Dynamics approach.In the parallel update-strategy, the agents act conceptually all at the same time in lock-step.This implies that they observe the same environment state during a time-step and actions of an agent are only visible in the next time-step -they are isolated from each other, see Figure 3.
Also, the authors [42] have shown the influence of different deterministic and non-deterministic elements in ABS on the dynamics and how the influence of non-determinism can completely break them down or result in different dynamics despite same initial conditions.This means that we want to rule out any potential source of non-determinism, which we achieve by keeping our implementation pure.This rules out the use of the IO Monad and thus any potential source of non-determinism under all circumstances because we would lose all compile time guarantees about reproducibility.Still we will make use of the Random Monad, which indeed allows side-effects but the crucial point here is that we restrict side-effects only to this type in a controlled way without allowing general unrestricted effects like in traditional object-oriented approaches in the field.
In the following, we derive a pure functional approach for an ABS of the SIR model in which we pose solutions to the previously mentioned problems.We start out with a first approach in Yampa and show its limitations.Then we generalise it to a more powerful approach, which utilises Monadic Stream Functions, a generalisation of FRP.Finally we add a structured environment, making the example more interesting and showing the real strength of ABS over other simulation methodologies like System Dynamics and Discrete Event Simulation1 .

Functional Reactive Programming
As described in the Section 2.2, Arrowized FRP [22] is a way to implement systems with continuous and discrete time-semantics where the central concept is the signal function, which can be understood as a process over time, mapping an input-to an outputsignal.Technically speaking, a signal function is a continuation which allows to capture state using closures and hides away the ∆t, which means that it is never exposed explicitly to the programmer, meaning it cannot be messed with.
As already pointed out, agents need to perceive time, which means that the concept of processes over time is an ideal match for our agents and our system as a whole, thus we will implement them and the whole system as signal functions.
4.1.1Implementation.We start by defining the SIR states as ADT and our agents as signal functions (SF) which receive the SIR Depending on the initial state we return the corresponding behaviour.Note that we are passing a random-number generator instead of running in the Random Monad because signal functions as implemented in Yampa are not capable of being monadic.
We see that the recovered agent ignores the random-number generator because a recovered agent does nothing, stays immune forever and can not get infected again in this model.Thus a recovered agent is a consuming state from which there is no escape, it simply acts as a sink which returns constantly Recovered: Next, we implement the behaviour of a susceptible agent.It makes contact on average with β other random agents.For every infected agent it gets into contact with, it becomes infected with a probability of γ .If an infection happens, it makes the transition to the Infected state.To make contact, it gets fed the states of all agents in the system from the previous time-step, so it can draw random contacts -this is one, very naive way of implementing the interactions between agents.
Thus a susceptible agent behaves as susceptible until it becomes infected.Upon infection an Event is returned, which results in switching into the infectedAgent SF, which causes the agent to behave as an infected agent from that moment on.When an infection event occurs we change the behaviour of an agent using the Yampa combinator switch, which is quite elegant and expressive as it makes the change of behaviour at the occurrence of an event explicit.Note that to make contact on average, we use Yampas occasionally function which requires us to carefully select the right ∆t for sampling the system as will be shown in results.
Note the use of iPre :: a → SF a a, which delays the input signal by one sample, taking an initial value for the output at time zero.The reason for it is that we need to delay the transition from susceptible to infected by one step due to the semantics of the switch combinator: whenever the switching event occurs, the signal function into which is switched will be run at the time of the event occurrence.This means that a susceptible agent could make a transition to recovered within one time-step, which we want to prevent, because the semantics should be that only one state-transition can happen per time-step.To deal with randomness in an FRP way, we implemented additional signal functions built on the noiseR function provided by Yampa.This is an example for the stream character and statefulness of a signal function as it allows to keep track of the changed randomnumber generator internally through the use of continuations and closures.Here we provide the implementation of randomBoolSF.drawRandomElemSF works similar but takes a list as input and returns a randomly chosen element from it: An infected agent recovers on average after δ time units.This is implemented by drawing the duration from an exponential distribution [5] with λ = 1 δ and making the transition to the Recovered state after this duration.Thus the infected agent behaves as infected until it recovers, on average after the illness duration, after which it behaves as a recovered agent by switching into recoveredAgent.As in the case of the susceptible agent, we use the occasionally function to generate the event when the agent recovers.Note that the infected agent ignores the states of the other agents as its behaviour is completely independent of them.What we need to implement next is a closed feedback-loop -the heart of every agent-based simulation.Fortunately, [9,29] discusses implementing this in Yampa.The function stepSimulation is an implementation of such a closed feedback-loop.It takes the current signal functions and states of all agents, runs them all in parallel and returns this step's new agent states.Note the use of notYet, which is required because in Yampa switching occurs immediately at t = 0.If we don't delay the switching at t = 0 until the next step, we would enter an infinite switching loop -notYet simply delays the first switching until the next time-step.Its first argument is the pairing-function, which pairs up the input to the signal functions -it has to preserve the structure of the signal function collection.The second argument is the collection of signal functions to run.The third argument is a signal function generating the switching event.The last argument is a function, which generates the continuation after the switching event has occurred.dpSwitch returns a new signal function, which runs all the signal functions in parallel and switches into the continuation when the switching event occurs.The d in dpSwitch stands for decoupled which guarantees that it delays the switching until the next time-step: the function into which we switch is only applied in the next step, which prevents an infinite loop if we switch into a recursive continuation.
Conceptually, dpSwitch allows us to recursively switch back into the stepSimulation with the continuations and new states of all the agents after they were run in parallel.

Results
. The dynamics generated by this step can be seen in Figure 4.
By following the FRP approach we assume a continuous flow of time, which means that we need to select a correct ∆t, otherwise we would end up with wrong dynamics.The selection of a correct ∆t depends in our case on occasionally in the susceptible behaviour, which randomly generates an event on average with contact rate following the exponential distribution.To arrive at the correct dynamics, this requires us to sample occasionally, and thus the whole system, with small enough ∆t which matches the frequency of events generated by contact rate.If we choose a too large ∆t, we loose events, which will result in wrong dynamics as can be seen in Figure 4a.This issue is known as under-sampling and is described in Figure 5.    5: A visual explanation of under-sampling and supersampling.The black dots represent the time-steps of the simulation.The red dots represent virtual events which occur at specific points in continuous time.In the case of undersampling, 3 events occur in between the two time steps but occasionally only captures the first one.By increasing the sampling frequency either through a smaller ∆t or supersampling all 3 events can be captured.
For tackling this issue we have two options.The first one is to use a smaller ∆t as can be seen in 4b, which results in the whole system being sampled more often, thus reducing performance.The other option is to implement super-sampling and apply it to occasionally, which would allow us to run the whole simulation with ∆t = 1.0 and only sample the occasionally function with a much higher frequency.

Discussion. We can conclude that our first step already introduced most of the fundamental concepts of ABS:
• Time -the simulation occurs over virtual time which is modelled explicitly, divided into fixed ∆t, where at each step all agents are executed.• Agents -we implement each agent as an individual, with the behaviour depending on its state.• Feedback -the output state of the agent in the current timestep t is the input state for the next time-step t + ∆t.• Environment -as environment we implicitly assume a fullyconnected network (complete graph) where every agent 'knows' every other agent, including itself and thus can make contact with all of them.• Stochasticity -it is an inherently stochastic simulation, which is indicated by the random-number generator and the usage of occasionally, randomBoolSF and drawRandomElemSF.• Deterministic -repeated runs with the same initial randomnumber generator result in same dynamics.This may not come as a surprise but in Haskell we can guarantee that property statically already at compile time because our simulation runs not in the IO Monad.This guarantees that no external, uncontrollable sources of non-determinism can interfere with the simulation.• Parallel, lock-step semantics -the simulation implements a parallel update-strategy where in each step the agents are run isolated in parallel and don't see the actions of the others until the next step.
Using FRP in the instance of Yampa results in a clear, expressive and robust implementation.State is implicitly encoded, depending on which signal function is active.By using explicit time-semantics with occasionally we can achieve extremely fine grained stochastics by sampling the system with small ∆t: we are treating it as a truly continuous time-driven agent-based system.
A very severe problem, hard to find with testing but detectable with in-depth validation analysis, is the fact that in the susceptible agent the same random-number generator is used in occasionally, drawRandomElemSF and randomBoolSF.This means that all three stochastic functions, which should be independent from each other, are inherently correlated.This is something one wants to prevent under all circumstances in a simulation, as it can invalidate the dynamics on a very subtle level, and indeed we have tested the influence of the correlation in this example and it has an impact.We left this severe bug in for explanatory reasons, as it shows an example where functional programming actually encourages very subtle bugs if one is not careful.A possible but not very elegant solution would be to simply split the initial random-number generator in sirAgent three times (using one of the splited generators for the next split) and pass three random-number generators to susceptible.A much more elegant solution would be to use the Random Monad which is not possible because Yampa is not monadic.
So far we have an acceptable implementation of an agent-based SIR approach.What we are lacking at the moment is a general treatment of an environment and an elegant solution to the random number correlation.In the next step we make the transition to Monadic Stream Functions as introduced in Dunai [34], which allows FRP within a monadic context and gives us a way for an elegant solution to the random number correlation.

Generalising to Monadic Stream Functions
A part of the library Dunai is BearRiver, a wrapper which reimplements Yampa on top of Dunai, which should allow us to easily replace Yampa with MSFs.This will enable us to run arbitrary monadic computations in a signal function, solving our problem of correlated random numbers through the use of the Random Monad.The question is now how to access this Random Monad functionality within the MSF context.For the function occasionally, there exists a monadic pendant occasionallyM which requires a MonadRandom type-class.Because we are now running within a MonadRandom instance we simply replace occasionally with occa-sionallyM.

Adding an environment
So far we have implicitly assumed a fully connected network amongst agents, where each agent can see and 'knows' every other agent.This is a valid environment and in accordance with the System Dynamics inspired implementation of the SIR model but does not show the real advantage of ABS to situate agents within arbitrary environments.Often, agents are situated within a discrete 2D environment [15] which is simply a finite N xM grid with either a Moore or von Neumann neighbourhood (Figure 6).Agents are either static or can move freely around with cells allowing either single or multiple occupants.
We can directly map the SIR model to a discrete 2D environment by placing the agents on a corresponding 2D grid with an unrestricted neighbourhood.The behaviour of the agents is the same but they select their interactions directly from the shared read-only environment, which will be passed to the agents as input.This allows agents to read the states of all their neighbours, which tells them if a neighbour is infected or not.To show the benefit over the System Dynamics approach and for purposes of a more interesting approach, we restrict the neighbourhood to Moore (Figure 6b).
We also implemented this spatial approach in Java using the well known ABS library RePast [30], to have a comparison with a state  of the art approach and came to the same results as shown in Figure 7.This supports, that our pure functional approach can produce such results as well and compares positively to the state of the art in the ABS field.
4.3.1 Implementation.We start by defining the discrete 2D environment for which we use an indexed two dimensional array.Each cell stores the agent state of the last time-step, thus we use the SIRState as type for our array data.Also, we re-define the agent signal function to take the structured environment SIREnv as input instead of the list of all agents as in our previous approach.As output we keep the SIRState, which is the state the agent is currently in.Also we run in the Random Monad as introduced before to avoid the random number correlation.Note that the environment is not returned as output because the agents do not directly manipulate the environment but only read from it.Again, this enforces the semantics of the parallel update-strategy through the types where the agents can only see the previous state of the environment and see the actions of other agents reflected in the environment only in the next step.
Note that we could have chosen to use a StateT transformer with the SIREnv as state, instead of passing it as input, with the agents then able to arbitrarily read/write, but this would have violated the semantics of our model because actions of agents would have become visible within the same time-step.
The implementation of the susceptible, infected and recovered agents are almost the same with only the neighbour querying now slightly different.
Stepping the simulation needs a new approach because in each step we need to collect the agent outputs and update the environment for the next next step.For this we implemented a separate MSF, which receives the coordinates for every agent to be able to update the state in the environment after the agent was run.Note that we need use mapM to run the agents because we are running now in the context of the Random Monad.This has the consequence that the agents are in fact run sequentially one after the other but because they cannot see the other agents actions nor observe changes in the shared read-only environment, it is conceptually a parallel update-strategy where agents run in lock-step, isolated from each other at conceptually the same time.--run agents sequentially but with shared, read-only environment ret <-mapM (`unMSF`env) sfs --construct new environment from all agent outputs for next step let (as, sfs') = unzip ret env' = foldr (\ (a, coord) envAcc -> updateCell coord a envAcc) env (zip as coords) sfsCoords' = zip sfs' coords cont = simulationStep sfsCoords' env' return (env', cont)) updateCell :: Disc2dCoord -> SIRState -> SIREnv -> SIREnv 4.3.2Results.We implemented rendering of the environments using the gloss library which allows us to cycle arbitrarily through the steps and inspect the spreading of the disease over time visually as seen in Figure 7.
Note that the dynamics of the spatial SIR simulation, which are seen in Figure 7b look quite different from the reference dynamics of Figure 2.This is due to a much more restricted neighbourhood which results in far fewer infected agents at a time and a lower number of recovered agents at the end of the epidemic, meaning that fewer agents got infected overall.

Discussion
. By introducing a structured environment with a Moore neighbourhood, we showed the ABS ability to place the heterogeneous agents in a generic environment, which is the fundamental advantage of an agent-based approach over other simulation methodologies and allows us to simulate much more realistic scenarios.
Note, that an environment is not restricted to be a discrete 2D grid and can be anything from a continuous N-dimensional space to a complex network -one only needs to change the type of the environment and agent input and provide corresponding neighbourhood querying functions.

Additional Steps
ABS involves a few more advanced concepts, which we don't fully explore in this paper due to lack of space.Instead we give a short overview and discuss them without presenting code or going into technical details.

Synchronous Agent
Interactions.Synchronous agent interactions are necessary when an arbitrary number of interactions between two agents need to happen instantaneously within the same time-step.The use-case for this are price negotiations between multiple agents where each pair of agents needs to come to an agreement in the same time-step [15].In object-oriented programming, the concept of synchronous communication between agents is implemented directly with method calls.We have implemented synchronous interactions in an additional step.We solved it pure functionally by running the signal functions of the transacting agent pair as often as their protocol requires but with ∆t = 0, which indicates the instantaneous character of these interactions.
4.4.2Event-Driven Approach.Our approach is inherently timedriven where the system is sampled with fixed ∆t.The other fundamental way to implement an ABS in general, is to follow an event-driven approach [28], which is based on the theory of Discrete Event Simulation [48].In such an approach the system is not sampled in fixed ∆t but advanced as events occur, where the system stays constant in between.Depending on the model, in an event-driven approach it may be more natural to express the requirements of the model.
In an additional step we have implemented a rudimentary eventdriven approach, which allows the scheduling of events.Using the flexibility of MSFs we added a State transformer to the monad stack, which allows queuing of events into a priority queue.The simulation is advanced by processing the next event at the top of the queue, which means running the MSF of the agent which receives the event.The simulation terminates if there are either no more events in the queue or after a given number of events, or if the simulation time has advanced to some limit.Having made the transition to MSFs, implementing this feature was quite straight forward, which shows the power and strength of the generalised approach to FRP using MSFs.

Conflicts in Environment.
The semantics of the agentbased SIR model allowed a straight-forward implementation of the parallel update-strategy.This is not easily possible when there could be conflicts in the environment e.g.moving agents where only a single one can occupy a cell.Most models in ABS [15] solve this by implementing a sequential update-strategy [42], where agents are run after another but can already observe the changes by agents run before them in the same time-step.To prevent the introduction of artefacts due to a specific ordering, these models shuffle the agents before running them in each step to average the probability for a specific agent to be run at a fixed position.
It is possible to implement a sequential update-strategy using the State Monad but functional programming might offer other conflict resolving mechanisms as well because of immutable data and its different nature of side-effects.One approach could be to still run the agents isolated from each other without a State Monad but in case of conflicts, to randomly select a winner and re-run other conflicting agents signal functions until there is no more conflict.As long as the underlying monadic context is robust to re-runs, e.g. the Random Monad, this is no problem.We argue that such an approach is conceptually and semantically cleaner and easier implemented in functional programming than in traditional object-oriented approaches.

RELATED WORK
The amount of research on using pure functional programming with Haskell in the field of ABS has been moderate so far.Most of the papers are related to the field of Multi Agent Systems and look into how agents can be specified using the belief-desire-intention paradigm [11,24,41].
The author of [4] investigated in his master thesis Haskells parallel and concurrency features to implement (amongst others) HLogo, a Haskell clone of the ABS simulation package NetLogo, where agents run within the IO Monad and make use of Software Transactional Memory for a limited form of agent-interactions.
A library for DES and SD in Haskell called Aivika 3 is described in the technical report [40].It is not pure, as it uses the IO Monad under the hood and comes only with very basic features for eventdriven ABS, which allows to specify simple state-based agents with timed transitions.
Using functional programming for DES was discussed in [24] where the authors explicitly mention the paradigm of FRP to be very suitable to DES.
A domain-specific language for developing functional reactive agent-based simulations was presented in [45].This language called FRABJOUS is human readable and easily understandable by domainexperts.It is not directly implemented in FRP/Haskell but is compiled to Yampa code which they claim is also readable.This supports that FRP is a suitable approach to implement ABS in Haskell.Unfortunately, the authors do not discuss their mapping of ABS to FRP on a technical level, which would be of most interest to functional programmers.
Object-oriented programming and simulation have a long history together as the former one emereged out of Simula 67 [10] which was created for simulation purposes.Simula 67 already supported Discrete Event Simulation and was highly influential for today's object-oriented languages.Although the language was important and influential, in our research we look into different approaches, orthogonal to the existing object-oriented concepts.
Lustre is a formally defined, declarative and synchronous dataflow programming language for programming reactive systems [17].
While it has solved some issues related to implementing ABS in Haskell it still lacks a few important features necessary for ABS.We don't see any way of implementing an environment in Lustre as we do in our approach in Section 4.3.Also the language seems not to come with stochastic functions, which are but the very building blocks of ABS.Finally, Lustre does only support static networks, which is clearly a drawback for ABS in general where agents can be created and terminated dynamically during simulation.
There exists some research [12,38,44] of using the functional programming language Erlang [3] to implement ABS.The language is inspired by the actor model [1] and was created in 1986 by Joe Armstrong for Eriksson for developing distributed high reliability software in telecommunications.The actor model can be seen as quite influential to the development of the concept of agents in ABS, which borrowed it from Multi Agent Systems [47].It emphasises message-passing concurrency with share-nothing semantics, which maps nicely to functional programming concepts.The mentioned papers investigate how the actor model can be used to close the conceptual gap between agent-specifications, which focus on messagepassing and their implementation.Further they also showed that using this kind of concurrency allows to overcome some problems of low level concurrent programming as well.Despite the natural mapping of ABS concepts to such an actor language, it leads to simulations, which despite same initial starting conditions, might result in different dynamics each time due to concurrency.

CONCLUSIONS
Our FRP based approach is different from traditional approaches in the ABS community.First it builds on the already quite powerful FRP paradigm.Second, due to our continuous time approach, it forces one to think properly of time-semantics of the model and how small ∆t should be.Third it requires one to think about agent interactions in a new way instead of being just method-calls.
Because no part of the simulation runs in the IO Monad and we do not use unsafePerformIO we can rule out a serious class of bugs caused by implicit data-dependencies and side-effects, which can occur in traditional imperative implementations.
Also we can statically guarantee the reproducibility of the simulation, which means that repeated runs with the same initial conditions are guaranteed to result in the same dynamics.Although we allow side-effects within agents, we restrict them to only the Random Monad in a controlled, deterministic way and never use the IO Monad, which guarantees the absence of non-deterministic side effects within the agents and other parts of the simulation.
Determinism is also ensured by fixing the ∆t and not making it dependent on the performance of e.g. a rendering-loop or other system-dependent sources of non-determinism as described by [35].Also by using FRP we gain all the benefits from it and can use research on testing, debugging and exploring FRP systems [33,35].
Also we showed how to implement the parallel update-strategy [42] in a way that the correct semantics are enforced and guaranteed already at compile time through the types.This is not possible in traditional imperative implementations and poses another unique benefit over the use of functional programming in ABS.

Issues
Currently, the performance of the system does not come close to imperative implementations.We compared the performance of our pure functional approach as presented in Section 4.3 to an implementation in Java using the ABS library RePast [30].We ran the simulation until t = 100 on a 51x51 (2,601 agents) with ∆t = 0.1 (unknown in RePast) and averaged 8 runs.The performance results make the lack of speed of our approach quite clear: the pure functional approach needs around 72.5 seconds whereas the Java RePast version just 10.8 seconds on our machine to arrive at t = 100.It must be mentioned, that RePast does implement an event-driven approach to ABS, which can be much more performant [28] than a time-driven one as ours, so the comparison is not completely valid.Still, we have already started investigating speeding up performance through the use of Software Transactional Memory [18,19], which is quite straight forward when using MSFs.It shows very good results but we have to leave the investigation and optimization of the performance aspect of our approach for further research as it is beyond the scope of this paper.
Despite the strengths and benefits we get by leveraging on FRP, there are errors that are not raised at compile time, e.g.we can still have infinite loops and run-time errors.This was for example investigated in [37] where the authors use dependent types to avoid some run-time errors in FRP.We suggest that one could go further and develop a domain specific type system for FRP that makes the FRP based ABS more predictable and that would support further mathematical analysis of its properties.Furthermore, moving to dependent types would pose a unique benefit over the traditional object-oriented approach and should allow us to express and guarantee even more properties at compile time.We leave this for further research.
In our pure functional approach, agent identity is not as clear as in traditional object-oriented programming, where there is a quite clear concept of object-identity through the encapsulation of data and methods.Signal functions don't offer this strong identity and one needs to build additional identity mechanisms on top e.g. when sending messages to specific agents.
We can conclude that the main difficulty of a pure functional approach evolves around the communication and interaction between agents, which is a direct consequence of the issue with agent identity.Agent interaction is straight-forward in object-oriented programming, where it is achieved using method-calls mutating the internal state of the agent, but that comes at the cost of a new class of bugs due to implicit data flow.In pure functional programming these data flows are explicit but our current approach of feeding back the states of all agents as inputs is not very general.We have added further mechanisms of agent interaction which we had to omit due to lack of space.

FURTHER RESEARCH
We see this paper as an intermediary and necessary step towards dependent types for which we first needed to understand the potential and limitations of a non-dependently typed pure functional approach in Haskell.Dependent types are extremely promising in functional programming as they allow us to express stronger guarantees about the correctness of programs and go as far as allowing to formulate programs and types as constructive proofs, which must be total by definition [2,27,43].
So far no research using dependent types in agent-based simulation exists at all.In our next paper we want to explore this for the first time and ask more specifically how we can add dependent types to our pure functional approach, which conceptual implications this has for ABS and what we gain from doing so.We plan on using Idris [6] as the language of choice as it is very close to Haskell with focus on real-world application and running programs as opposed to other languages with dependent types e.g.Agda and Coq.
We hypothesize that dependent types could help ruling out even more classes of bugs at compile time and encode invariants and model specifications on the type level, which implies that we don't need to test them using e.g.property-testing with QuickCheck.This would allow the ABS community to reason about a model directly in code.We think that a promising approach is to follow the work of [7,8,16] in which the authors utilize GADTs to implement an indexed monad, which allows to implementation correct-byconstruction software.
• In the SIR implementation one could make wrong statetransitions e.g. when an infected agent should recover, nothing prevents one from making the transition back to susceptible.
Using dependent types it should be possible to encode invariants and state-machines on the type level, which can prevent such invalid transitions already at compile time.This would be a huge benefit for ABS because many agent-based models define their agents in terms of state-machines.• An infected agent recovers after a given time -the transition of infected to recovered is a timed transition.Nothing prevents us from never doing the transition at all.With dependent types we should be able to encode the passing of time in the types and guarantee on a type level that an infected agent has to recover after a finite number of time steps.• In more sophisticated models agents interact in more complex ways with each other e.g. through message exchange using agent IDs to identify target agents.The existence of an agent is not guaranteed and depends on the simulation time because agents can be created or terminated at any point during simulation.Dependent types could be used to implement agent IDs as a proof that an agent with the given id exists at the current time-step.This also implies that such a proof cannot be used in the future, which is prevented by the type system as it is not safe to assume that the agent will still exist in the next step.• In our implementation, we terminate the SIR model always after a fixed number of time-steps.We can informally reason that restricting the simulation to a fixed number of timesteps is not necessary because the SIR model has to reach a steady state after a finite number of steps.This means that at that point the dynamics won't change any more, thus one can safely terminate the simulation.Informally speaking, the reason for that is that eventually the system will run out of infected agents, which are the drivers of the dynamic.
We know that all infected agents will recover after a finite number of time-steps and that there is only a finite source for infected agents which is monotonously decreasing.Using dependent types it might be possible to encode this in the types, resulting in a total simulation, creating a correspondence between the equilibrium of a simulation and the totality of its implementation.Of course this is only possible for models in which we know about their equilibria a priori or in which we can reason somehow that an equilibrium exists.
fallingMass :: Double -> Double -> SF () Double fallingMass p0 v0 = proc _ -> do v <-arr (+v0) <<< integral -< (-9.8) p <-arr (+p0) <<< integral -< v returnA -< p newtype MSF m a b = MSF {unMSF :: MSF m a b -> a -> m (b, MSF m a b)} MSFs are also arrows, which means we can apply arrowized programming with Patersons do-notation as well.MSFs are implemented in Dunai, which is available on Hackage.Dunai allows us to apply monadic transformations to every sample by means of combinators like arrM :: Monad m ⇒ (a → m b) → MSF m a b and arrM_ :: Monad m ⇒ m b → MSF m a b.A part of the library Dunai is BearRiver, a wrapper, which re-implements Yampa on top of Dunai, which enables one to run arbitrary monadic computations in a signal function.BearRiver simply adds a monadic parameter m to each SF, which indicates the monadic context this signal function runs in.

Figure 1 :
Figure 1: States and transitions in the SIR compartment model.

Figure 3 :
Figure 3: Parallel, lock-step execution of the agents.

Figure 4 :
Figure 4: FRP simulation of agent-based SIR showing the influence of different ∆t.Population size of 1,000 with contact rate β = 1 5 , infection probability γ = 0.05, illness duration δ = 15 with initially 1 infected agent.Simulation run for 150 time-steps with respective ∆t.

Figure
Figure5: A visual explanation of under-sampling and supersampling.The black dots represent the time-steps of the simulation.The red dots represent virtual events which occur at specific points in continuous time.In the case of undersampling, 3 events occur in between the two time steps but occasionally only captures the first one.By increasing the sampling frequency either through a smaller ∆t or supersampling all 3 events can be captured.

4. 2 . 1
Identity Monad.We start by making the transition to Bear-River by simply replacing Yampas signal function by BearRivers', which is the same but takes an additional type parameter m, indicating the monadic context.If we replace this type-parameter with the Identity Monad, we should be able to keep the code exactly the same, because BearRiver re-implements all necessary functions we are using from Yampa.We simply re-define the agent signal function, introducing the monad stack our SIR implementation runs in: type SIRMonad = Identity type SIRAgent = SF SIRMonad [SIRState] SIRState 4.2.2Random Monad.Using the Identity Monad does not gain us anything but it is a first step towards a more general solution.Our next step is to replace the Identity Monad by the Random Monad, which will allow us to run the whole simulation within the Random Monad with the full features of FRP, finally solving the problem of correlated random numbers in an elegant way.We start by re-defining the SIRMonad and SIRAgent: type SIRMonad g = Rand g type SIRAgent g = SF (SIRMonad g) [SIRState] SIRState occasionallyM :: MonadRandom m => Time -> b -> SF m a (Event b) --can be used through the use of arrM and lift randomBoolM :: RandomGen g => Double -> Rand g Bool --this can be used directly as a SF with the arrow notation drawRandomElemSF :: MonadRandom m => SF m [a] a 4.2.3Discussion.Running in the Random Monad solved the problem of correlated random numbers and elegantly guarantees us that we won't have correlated stochastics as discussed in the previous section.In the next step we introduce the concept of an explicit discrete 2D environment.

Figure 6 :
Figure 6: Common neighbourhoods in discrete 2D environments of Agent-Based Simulation.

( a )Figure 7 :
Figure 7: Simulating the agent-based SIR model on a 21x21 2D grid with Moore neighbourhood (Figure 6b), a single infected agent at the center and same SIR parameters as in Figure 2. Simulation run until t = 200 with fixed ∆t = 0.01.Last infected agent recovers around t = 194.The susceptible agents are rendered as blue hollow circles for better contrast.
• States -agents encapsulate some state, which can be accessed and changed during the simulation.