BharatSim Documentation

Welcome to the documentation for BharatSim!

BharatSim is a collaborative project between Ashoka University and Thoughtworks, funded by the Bill & Melinda Gates Foundation. Its ongoing development at Ashoka University is funded by the Mphasis F1 Foundation.

BharatSim was originally designed to run decision-critical scenarios for India during the COVID-19 pandemic. Standard compartmental models that are widely used to model the spread of infectious diseases largely ignore individual level granularity, since they mainly consider a well-mixed population. Such models are ill-suited for predicting and analysing complex real-world phenomena. Realistic systems involve interactions between individuals with different personal attributes (age, weight, height, and comorbidities) and geographies. These interactions can potentially lead to emergent phenomena, especially in scenarios such as pandemics where different individuals are affected differently, based on their attributes.

An agent-based approach allow modellers to account for individual-level heterogeneity, and can thus be used to simulate intricate scenarios of varying complexity in which agents interact with each other and their environment. The results of such simulations could, in principle, be used to guide policy-level interventions like lockdowns and vaccination drives.

BharatSim’s vision is to build a simulation framework that is distributed, multi-scale, and agent-based for use by the scientific community. It aims to provide the user with a synthetic population tailored to the Indian context, and a modelling framework that can easily be customised to address a range of questions. The framework has two components:

  • Simulation engine: Given a Synthetic Population, the simulation engine can support large-scale simulations with multi-million agents, incorporating daily behaviours and policy-level interventions. It structures this information by treating agents and locations as nodes, and the relations between them as edges on a graph network. This way modellers can analyse how the network structure of the population affects the spread of the disease. The simulation engine is written in Scala, and allows modellers to directly specify their models using a domain-specific, high-level language. This domain-specific language is also based on the Scala programming language, thus allowing modellers to extend their knowledge of Scala when creating their models.

  • Visualization engine: BharatSim also includes a dashboard visualisation tool that allows the user to create and view multiple types of graphs at the same time. The visualisation engine was developed to plot and analyse the outputs generated by the simulation engine, but can be used as an independent application as well. It possesses a customisable dashboard system that is flexible, allowing users to produce graphs that study both the temporal and spatial variation of disease spread.

In this tutorial documentation, we look be looking at BharatSim in detail, introducing the novice programmer to both components. The tutorial assumes a familiarity with (or at least an eagerness to learn) Object-oriented Programming. Prior familiarity with Scala or Java is desirable, but not a prerequisite.

Introduction

Simulating hypothetical scenarios is a very useful way to analyse and predict the behaviour of complex phenomena such as the spread of COVID-19 through a network of individuals. An infectious disease outbreak is heavily influenced by factors pertaining to both the disease (like the transmissibility of the disease, or how it is affected by vaccines), and the network on which it spreads (such as the age structure, population density, network contacts, and geography). These factors can be responsible for heterogeneity in the size and severity of outbreaks in different locations that might otherwise appear similar.

For this reason the results from models developed for other countries may not be able to achieve the desirable level of predictive ability. Such models might not account for factors pertinent to India – its unique demographics, the state (or absence) of health care, stratified social structures, and complex ecological gradients. The inputs for such factors must be collated from multiple sources, including imperfect and poorly usable government data, which requires both local knowledge and considerable experience in data analysis and interpretation. The result is a pressing need for an epidemic modelling framework tailored to Indian needs, with potential applications beyond its immediate use for COVID-19 modelling.

Given these design requirements, building and implementing detailed agent-based models is an unsurprisingly daunting task. A useful model should ideally predict impending potential surges in healthcare requirements, allowing for timely resource allocation. It should also allow the modeler to compare the merits and consequences of different non-pharmaceutical interventions at different stages of the epidemic. As a specific example, a useful model might help the modeler untangle how asymptomatic infections play a role in the spread of the disease.

BharatSim has been developed with 3 goals in mind:

  1. Flexibility: Since it is a social simulation framework, a researcher can develop models for a variety of disciplines including epidemiology, economics, climate science.

  2. Scalability: It can keep track of millions of agents in an efficient manner.

  3. Customisation: It has been developed to suit India’s needs, via its support for synthetic population.

BharatSim has been designed keeping these design requirements and challenges of the Indian context in mind. BharatSim’s vision is to build an India scale agent-based framework that would enable modellers from different disciplines – ranging from epidemiology, disaster management, and economics – to advise policy makers and decision makers across institutions.

Synthetic Population

A synthetic population is a simplified individual-level representation of the actual population. This means that while every person is represented individually in it, not all of their attributes are included (for example, hair colour or shoe-size are deemed to be irrelevant for modelling epidemic spread, and are thus ignored, while the presence of comorbidities like diabetes would be included). As such, a synthetic population does not aim to be identical to the actual population, but instead attempts to match its various statistical distributions and correlations, thereby being sufficiently close to the true population to be used in modelling.

In the table below, you can see an example of a section of a synthetic population. Each row represents an individual with a unique ID, as well as certain attributes. These attributes could be related to the individual themselves (like their gender, age, and height and so on), or their network (details pertaining to their homes, workplaces, and possibly schools). Additionally, the population could also contain information regarding the individual’s comorbidities (for example, whether they have diabetes or other preexisting conditions), if this is deemed relevant to the modelling exercise.

sample_synthetic_population.csv

Agent_ID

SexLabel

Age

Height

Weight

JobLabel

StateLabel

District

AdminUnitName

AdminUnitLatitude

AdminUnitLongitude

HHID

H_Lat

H_Lon

WorkPlaceID

W_Lat

W_Lon

school_id

school_lat

school_long

public_place_id

public_place_lat

public_place_long

essential_worker

Adherence_to_Intervention

M_High_BP

M_Diabetes

99802077568

Female

76

147.37

49.05

Construction

Maharashtra

Pune

Mohammadwadi-Kausar Baug

18.47477

73.89257

99801684702

18.46709

73.90603

2001000003467

18.4679

73.89859

0

3001000000062

18.45035

73.87068

0

1

0

0

99801380107

Male

37

162.03

57.94

Ag labour

Maharashtra

Pune

Nagpur Chawl-Phule Nagar

18.55893

73.87609

99801248473

18.55952

73.87877

2001000006630

18.58283

73.91661

0

3001000001044

18.52699

73.83451

1

0

0

0

99802408169

Male

6

111.21

23.13

Student

Maharashtra

Pune

Kharadi-Chandan Nagar

18.56355

73.93256

99800525921

18.54846

73.94971

0

2001000002070

18.55683

73.94757

3001000000274

18.54904

73.9491

0

1

0

0

99800402683

Female

37

152.65

52.61

Ag labour

Maharashtra

Pune

Karve Nagar

18.49149

73.82173

99800473441

18.48539

73.82129

2001000006876

18.53875

73.92594

0

3001000000650

18.48382

73.79731

1

0

0

0

99800824202

Female

35

150.92

52.42

Homebound

Maharashtra

Pune

Deccan Gymkhana-Model Colony

18.51845

73.83391

99800895513

18.52335

73.85339

0

0

3001000000236

18.54026

73.91186

0

0

0

0

99801178045

Female

50

151.51

50.1

Ag labour

Maharashtra

Pune

Shanivar Peth-Sadashiv Peth

18.51944

73.85194

99801142021

18.51388

73.84935

2001000000636

18.50785

73.84921

0

3001000001403

18.51024

73.84731

0

0.9

0

0

99802883562

Female

14

139.25

37.81

Student

Maharashtra

Pune

Viman Nagar-Somnath Nagar

18.55818

73.92165

99801268501

18.57398

73.9285

0

2001000000159

18.52152

73.92834

3001000000633

18.56121

73.93799

0

0.8

0

0

99800390824

Male

19

157.71

50.84

Plantation lab

Maharashtra

Pune

Salisbury Park-Maharshi Nagar

18.49162

73.8661

99800461238

18.49814

73.87331

2001000011846

18.49934

73.8511

0

3001000000528

18.46247

73.85122

0

0.4

0

0

99800144374

Female

50

152.6

50.58

Homebound

Maharashtra

Pune

Dhanakvadi-Ambegaon Pathar

18.45878

73.84264

99800193269

18.46661

73.8561

0

0

3001000001414

18.46532

73.89147

0

0.9

0

0

99800944765

Female

45

151.99

50.99

Homebound

Maharashtra

Pune

Kharadi-Chandan Nagar

18.56355

73.93256

99800978186

18.54808

73.94093

0

0

3001000001402

18.5684

73.82062

0

0.9

0

0

All of these attributes are strongly correlated with each other and a good synthetic population will ideally be able reproduce the correlations that occur in the real world. However, this is a monumental task; real world data is complex, and often contains many artifacts that need to be addressed.

A Basic Introduction to Epidemic Modelling

A standard technique to model the spread of epidemics is through a compartmental model. In such models, we divide the population into different categories or compartments according to the infection stage (such as “Susceptible”, “Infected”, and “Recovered”, for example). Since individuals can proceed through the different stages of the infection over time, we also incorporate transition rates between these compartments. In the following section we offer an introduction to a basic compartmental models, starting with the simple but surprisingly effective SIR model.

The SIR Model

The disease progression in the SIR Model

Compartments of the SIR model

The individuals in the SIR model can be in one of the three compartments at any given time, (S)usceptible, (I)nfected, or (R)ecovered. Infected individuals can “infect” the Susceptible, causing them to transition to the Infected compartment, from which they will eventually recover. In what follows we will use the letters S, I, and R to represent both the individual compartments as well as the numbers in each of these compartments. The context should make it clear what we are referring to. In addition, we make some assumptions:

  1. Susceptible individuals can only leave this compartment if they have been infected, and the total number of individuals in the population is conserved. In other words, we ignore the effects of birth and death rates.

  2. The rate of change of \(S(t)\), the number of susceptibles, depends on the number of individuals currently susceptible, the number currently infected, and a parameter that governs the amount of contact between susceptibles and infected.

  3. We also assume that some fixed fraction (say, \(\gamma\)) of the infected group will recover on any given day. Keep in mind that when we say someone has “Recovered”, what we mean (in this model) is that they are no longer infectious or susceptible to the disease and therefore cannot contribute to its spread.

Mean-field or “Well-mixed” solution to this model

This system is modelled by a set of coupled, nonlinear, first order ordinary differential equations. S, I, R are the variables, while \(\beta\) and \(\gamma\) are the parameters of this model. The terms on the left hand side represent rates of change. In this case, change in population of S,I,R per day.

\[\begin{split} \begin{equation} \begin{aligned} \dv{S}{t} &= - \beta\, S(t)\, I(t)\\[10pt] \dv{I}{t} &= + \beta \, S(t) \, I(t) - \gamma \, I(t)\\[10pt] \dv{R}{t} &= + \gamma \, I(t) \end{aligned} \tag{1} \label{1} \end{equation} \end{split}\]

Exercise

  1. Why is there a negative sign in the equation for \(S(t)\)?

  2. Convince yourself that \(\beta\) represents the average number of contacts per unit time between an infected individual and a susceptible individual.

  3. Show that the total population \(N = S+I+R\) is fixed (i.e. that this is a closed system). This is often a very good check to see if your simulation is running correctly!

By inspecting the equations, you can see that the parameters \(\beta\) and \(\gamma\) differ in nature. \(\gamma\) is just the fraction of infected individuals recovering per day. \(\beta\), on the other hand, is a way of quantifying assumption no. 2 from earlier. We want to express that the probability of infection depends on the contact between the infected and the susceptible. This amount of contact is determined by multiplying S and I. We call \(\beta\) the transmission coefficient.

\[\text{Transmission coefficient}\beta = \text{Contact probability } \times \text{ Transmission rate},\]

Each infected individual meets a fraction of the total number of susceptibles and infects them at some “transmission rate”. The “contact probability” decides how the transmission scales with the population. Two specific types of scaling can be chosen, depending on the type of disease.

  1. Contact probability = 1. In this case, we assume that infected individuals meet all the susceptibles in one day. Therefore, the larger the population, the faster the spread of disease. This can be used to model a large number of people crowded in an enclosed space. Such a scaling is often appropriate in the case of plant and animal diseases. This is known as density dependent scaling.

  2. Contact probability = \(\mathbf{1/N}\). In this case, we assume the infected individual meets a constant number of people (say \(\lambda_S\)) every day. Of these, the fraction of susceptibles is \(S/N\). This assumption holds for most human diseases, where contact is determined by social factors. This is known as frequency dependent scaling.

The distinction between these two types of scaling only occurs when the total population size is not fixed, since otherwise the factor of \(N\) could simply be absorbed into the definition of \(\beta\). In what follows, we will be using frequency dependent scaling. As a result, Equation (\ref{1}) above can be rewritten as:

\[\begin{split}\begin{equation} \begin{aligned} \dv{S}{t} &= - \frac{\lambda_S}{N} S I\\[10pt] \dv{I}{t} &= + \frac{\lambda_S}{N} S I - \lambda_I I\\[10pt] \dv{R}{t} &= + \lambda_I I \end{aligned} \tag{2} \label{2} \end{equation}\end{split}\]

Exercise

  1. Show that, because of the extra factor of \(N\), the solutions to this differential equation will be independent of \(N\). (Hint: Rewrite the equations in terms of \(s(t) = S(t)/N\), \(i(t) = I(t)/N\), and \(r(t) = R(t)/N\).)

  2. It turns out that finding the exact solution to this equation as a function of time is not possible. However, an implicit solution can be found. By dividing the first and third equations, show that

    \[s(t) = s(0)\,e^{-R_0 r(t)},\]

    where the quantity \(R_0\) is called the reproductive ratio. Show that \(R_0\) is independent of \(N\).

  3. Now, repeat the previous question for the system denoted by Equation (\ref{1}). Show that in that case, \(R_0\) depends on the total population \(N\).

The above differential equations represented in Equation (\ref{2}) can be easily solved numerically to obtain the curves for \(S\), \(I\), and \(R\) shown in the figure below.

The deterministic solution to the differential equations of the SIR Model

Of course, these solutions are deterministic. This is because we assume the transition rates between the compartments are fixed, the population is well-mixed and we treat all individuals as identical.

Note

A well-mixed population is one in which any infected individual has a probability of contacting any susceptible individual that can be approximated reasonably well by the average probability of susceptible-infected interaction. This is often the most problematic assumption, but is easily relaxed in more complex models.

Stochastic solutions to this model

What happens when the population doesn’t behave in a “well-mixed” manner? For example, consider a population where individuals move between their homes and work-places. In this case, all individuals might not have the same number of contacts. Some individuals might work in high-density workplaces and come in contact with many more individuals and spread the disease at a faster rate than others. The well-mixed scenario also assumes that everyone in a population is identical. We might want to account for the heterogeneity of individuals in the population: some agents might intrinsically be more likely to get infected than others. And lastly we might also want to implement different interventions like a lockdown where only certain agents are allowed to move, and not others.

For all of the above cases, the well-mixed system is inadequate since it assumes that all individuals are identical and indistinguishable. To get around these limitations, one approach is to treat each individual as a separate agent with attributes. These heterogeneous agents interact with each other, spreading the infection. However, keeping track of individual agents is computationally very resource-intensive, even if the questions we can answer are broader.

One of the most well-known methods to implement such simulations is the Gillespie Algorithm. Our framework uses a much simpler discrete time approximation of this method. (The steps for this algorithm are outlined in the box below.) We first consider that all the individuals are in a single location, i.e. everyone is in contact with everyone else. However, in the next section, we will relax this assumption and allow for networks of individuals to be formed. The basic idea is as follows:

Algorithm

  1. Divide the total time into steps of $\Delta t$, and at every time-step we loop over all agents.

  2. If the agent is susceptible, we compute the number of infected individuals who could potentially infect them (”$I$”). Then, with some probability $$p_\text{SI} = \lambda_S\frac{I}{N}\Delta t,$$ we transition them to the infected compartment.

  3. If the agent is already infected, we transition them to the recovered compartment with a probability $$p_\text{IR} = \lambda_I\,\Delta t.$$

  4. If they have recovered, do nothing.

  5. Repeat the entire process until there are no more infected individuals, or the total time has elapsed.

The results of such a stochastic simulation are shown in the figure below. Each faintly visible curve represents a realisation of the stochastic algorithm starting from the same initial conditions. As you can see, the progress of the disease is no longer deterministic. However (in case all agents are in a single location), the average over all of these stochastic runs results in the “well-mixed” solutions (boldly visible curves). Comparing the average of several stochastic runs to the deterministic solution is one way to check your code.

The stochastic solutions to the differential equations of the SIR Model

Introducing a network structure

So far, we have considered a population in which every agent is connected with every other agent. However, in real populations, individuals are generally only in contact with a small section of the total population. These contact structures depend heavily on the locations the individual frequents, and other agents who frequent this location, and can – in principle – depend on other attributes like their age, gender, and socio-economic class. In order to implement this heterogeneity, we will now need to implement an underlying network structure to our model, with individuals being associated with certain locations, and interacting with other individuals who are associated with those locations.

Consider the semi-realistic scenario where agents can move between different locations like “work” and “home”. We can build a network that accounts for the relationship between agents and locations over time.

A schematic of a model network

By adding this complexity, we can see that the number of individuals at any given location is not fixed. So we modify our algorithm such that we treat each location as a well-mixed case looping over the instantaneous population at the given location. The modified algorithm is shown in the box below.

Algorithm

  1. Divide the total time into steps of \(\Delta t\), and at every time-step we loop over all agents.

  2. If the agent is susceptible, we compute the number of infected individuals who could potentially infect them (”\(I\)”) in their current location. Then, with some probability

    \[p_\text{SI} = \lambda_S\frac{I}{N}\Delta t,\]
    we transition them to the infected compartment.

  3. If the agent is already infected, we transition them to the recovered compartment with a probability

    \[p_\text{IR} = \lambda_I\,\Delta t.\]

  4. If they have recovered, do nothing.

  5. Repeat the entire process until there are no more infected individuals, or the total time has elapsed.

Sojourn times in different compartments

How long does an individual spend in an infected state or compartment? This duration, also known as the sojourn time, isn’t a fixed value, so it’s important to know the distribution of these sojourn times. In Markov-Chain-Monte-Carlo simulations like ours, the sojourn times are exponentially distributed. First we will explain why the distribution is exponential, which can allow us to change our algorithm for non-exponential cases.

We discretise time into small intervals or time-steps of \(\Delta t\), so that the total time \(T = N \Delta t\), where \(N\) is the total number of time-steps. In every small step \(\Delta t\) an individual in our disease model has a small probability, \(p\), of exiting the compartment. The value of \(p\) depends on factors such as the number of individuals in that compartment at that time, and so on.

To simulate \(p\), we draw a random number \(r\) from a uniform distribution. If \(r < p\), the individual exits the compartment. This is equivalent to saying:

\[\begin{equation}\texttt{The probability of exiting the compartment is p}.\end{equation}\]

This is the only role of the uniform distribution here.

Now, we ask ourself the following question: what is the probability that an individual will leave a compartment at some time \(t\)? In our scheme, this is equivalent to asking what the probability is of the event occurring between time \(t\) and \(t + \Delta t\). We will call this probability \(P(t)\Delta t\). We’ve also discretized time into units of \(\Delta t\), so we can define \(t = n\Delta t\), where \(n\) is the number of steps of \(\Delta t\).

This allows us to reframe the problem as this: what is the probability that the individual will leave the compartment after exactly \(n\) steps? This is given by:

\[\text{Prob. that event occurs exactly after $n$ steps} = (1-p)^n p.\]

Using the fact that \(p = \lambda \Delta t\), and that \(\Delta t = T/N\), we can show that:

\[\begin{split} \begin{align} \mathcal{P}_{\Delta t}(t)\dd t &= (1 - \lambda \Delta t)^{t/\Delta t} \lambda \Delta t\\ &= \left(1 - \frac{\lambda T}{N} \right)^{t N/T} \lambda \Delta t \end{align} \end{split}\]

To find the true probability density, we need to take the limit \(N\to\infty\), i.e. \(N/T \to \infty\), and so

\[\mathcal{P}(t)\dd t = \lim_{N/T \to\infty} \left(1 - \frac{\lambda T}{N} \right)^{t N/T} \lambda \Delta t = \lambda e^{-\lambda t} \dd t, \]

so that the probability distribution is:

\[\mathcal{P}(t) = \lambda e^{-\lambda t}.\]

Exercise

  1. Show that the probability that an individual exits the compartment $I$ exactly after some time $t = n \Delta t$ is given by: $$\mathcal{P}(t)\,\Delta t = p_\text{IR} (1-p_\text{IR})^n.$$

  2. Using the fact that $p=\lambda_I\,\Delta t$, and $\Delta t = T/N_\text{steps}$ (Where $T$ is the total simulation time and $N_\text{steps}$ the total number of steps): $$\mathcal{P}(t)\,\Delta t = \left(1 - \frac{\lambda_I T}{N_\text{steps}}\right)^{t N_\text{steps}/T}\, \lambda_I \, \Delta t.$$

  3. Next, take the limit $N_\text{steps}\to\infty$, and $\Delta t\to 0$, keeping $N_\text{steps}\,\Delta t = T$, and argue that $$\mathcal{P}(t) = \lambda_I \, e^{-\lambda_I t}.$$ The sojourn times in the infected compartment are exponentially distributed, with mean $\tau_I = 1/\lambda_I$! In other words, this is a succession of Poisson processes – characteristic of Markov Chain Monte Carlo processes: happens every time the probability of a transition is independent of the history.

  4. Can the same argument be used to say that the sojourn time in the susceptible compartment is $\tau_S = 1/\lambda_S$? Explain.

Other Compartmental Models

The SEIR Model

The SEIR model is a simple extension of the basic SIR model to include the Exposed state to account for disease pathologies that go through a latent phase where an individual hasn’t started infecting others despite being exposed to the disease. The Exposed compartment (E) represents this incubation period for the disease. The infected individuals expose the Susceptible individuals (S) to the disease, moving the latter into the Exposed (E) compartment before they are moved to the Infected (I) compartment. From the infected compartment they will be Removed (R) eventually. Here we use the “Removed” compartment to indicate that individuals are removed from contributing to the disease spread. This could mean permanently recovered, immune, or dead. The diagram below shows how the individuals move through each compartment in this model.

The disease progression in the SEIR Model

The rate at which the disease is transmitted from an infected agent to a susceptible agent is represented by \({\lambda_S}\) (transmission rate). The incubation rate, \({\lambda_E}\), is the rate of Exposed individuals becoming infectious. The average time an individual spends in the Exposed compartment is given by \({1/\lambda_E}\). Lastly, \({\lambda_I}\) represents the rate of removal of infected individuals from the Infected compartment.

In a closed population with no births or deaths, the SEIR model can be defined using a set of coupled non-linear differential equations described below:

\[\begin{split}\begin{aligned} \dv{S}{t} &= -\lambda_S \frac{SI}{N} \\[10pt] \dv{E}{t} &= \lambda_S \frac{SI}{N} - \lambda_E E \\[10pt] \dv{I}{t} &= \lambda_E E - \lambda_I I \\[10pt] \dv{R}{t} &= \lambda_I I \end{aligned}\end{split}\]

where the total population,

\[N = S + E + I + R\]

Introducing the incubation period does not change the total number of infections, but prolongs the duration of the epidemic. The graphs below show simple SEIR models with incubation periods 5 and 10 days respectively.

_images/seir2.png _images/seir.png

The above equations can be solved numerically to get deterministic results but, as explained in The SIR Model, we can also solve it stochastically using a similar algorithm.

In the algorithm, if the agent is Susceptible, we compute the number of infected individuals they come in contact with who could potentially infect them (\(I\)). Then, during each time step \({\Delta t}\), they are transferred to the Exposed compartment, with some probability,

\[P_\text{SE} = \lambda_S \frac{I}{N}\Delta t\]

Individuals from the Exposed compartment are transferred to the Infected compartment with the probability,

\[P_\text{EI} = \lambda_E \Delta t\]

If the agent is already infected, we transition them to the Removed compartment with a probability

\[P_\text{IR} = \lambda_I \Delta t.\]

Similarly to the SIR model, the average of the stochastic solutions come close to the mean field solution, as seen in the figure below.

The stochastic solutions to the differential equations of the SEIR Model

\(\lambda_S=3/35 \text{ day}^{-1}\), \(\lambda_E=1/14 \text{ day}^{-1}\), \(\lambda_I=1/28 \text{ day}^{-1}\), \(N=10,000\).

Initial conditions: 1% Infected, 99% Susceptible

The SAIR Model

In the models we have described so far, we have not distinguished between the types of infected individuals in the population. In many real-world situations, however, we might need to make this distinction. For example, the disease progression of mildly infected individuals might be very different from severely infected individuals. Additionally, the existence of asymptomatic individuals may be crucial in studying the efficiency of non-pharmaceutical interventions such as quarantining individuals.

Let us begin by examining a very basic case: a generalisation of the The SIR Model to include both symptomatic (\(I\)) and asymptomatic (\(A\)) individuals. Such a distinction might be important to study the spread of epidemics like COVID-19, especially because asymptomatic individuals are much more likely to spread the disease as they are hard to indentify without extensive testing and contact tracing. From these compartments the individuals move to the Removed (\(R\)) compartment, at rates \(\lambda_A\) and \(\lambda_I\) respectively, as shown in the disease progression below.

In the models described so far, we have assumed that all infections are equal in their severity or intensity. However, in many real-world diseases, the disease manifests itself differently across individuals. For example, some individuals might be mildly infected, some might mount a severe symptomatic response, while others might be entirely asymptomatic. Accounting for these different types of infections is crucial to accurately model the disease progression through the population. Let us begin by examining a very basic case: a generalisation of the The SIR Model to include both symptomatic (\(I\)) and asymptomatic (\(A\)) individuals. Such a distinction might be important to study the spread of epidemics like COVID-19 and guide understanding and interventions. For instance, asymptomatic individuals might end up spreading the disease far and wide since they are hard to identify without extensive testing and contact tracing. Similarly, if we can predict the trajectory of severe cases, healthcare measures can be taken appropriately. From both these compartments the individuals move to the Removed (\(R\)) compartment, at rates \(\lambda_A\) and \(\lambda_I\) respectively, as shown in the disease progression below.

The disease progression in the SAIR Model

Initially, we can assume that both asymptomatic and symptomatic individuals infect susceptibles with equal capacity. (In reality, this capacity depends on various nuanced aspects of the disease pathology and the network of interacting individuals. However, we won’t delve into those details here.) We call this transition rate out of \(S\), \(\lambda_S\), as before.

However, now a branching event can occur. Once infected, a susceptible person could either move to \(A\) or \(I\). We thus define another quantity, \(\gamma\), which is the fraction of the infected individuals who are asymptomatic. The individuals then transit out of \(A\) or \(I\) with rates \(\lambda_A\) or \(\lambda_I\) respectively. The set of coupled non-linear differential equations that defines the SAIR model in a closed population are:

\[\begin{split}\begin{aligned} \dv{S}{t} &= -\frac{\lambda_S}{N} S\left(A + I\right) \\[10pt] \dv{A}{t} &= \gamma \frac{\lambda_S}{N} S \left(A + I\right) - \lambda_A A \\[10pt] \dv{I}{t} &= (1-\gamma) \frac{\lambda_S}{N} \left(A+I\right) - \lambda_I I \\[10pt] \dv{R}{t} &= \lambda_A A+ \lambda_I I \end{aligned}\end{split}\]

where, just as before, the total population is constant:

\[N = S + A + I + R.\]

Exercise

Convince yourself that if \(\lambda_A = \lambda_I\), this model effectively reduces to a simple \(SIR\) model. In this case the distinction between the asymptomatics and symptomatics is merely cosmetic.

Sample run for the SAIR Model

Modelling the transitions in the SAIR model is a little bit more involved than in the SIR model, though the basic principle is the same.

Warning

One might naively imagine that we could simply write:

\[\begin{split}P_\text{SA} &= \gamma \lambda_{S} \left(\frac{A+I}{N}\right) \Delta t,\\ P_\text{SI} &= (1-\gamma) \lambda_{S} \left(\frac{A+I}{N}\right) \Delta t,\end{split}\]

and draw two random numbers \(r_1\) and \(r_2\) to check if \(P_\text{SA}\) or \(P_\text{SI}\) occur. However, this is not strictly correct. The transitions from \(S\) to \(A\) and from \(S\) to \(I\) are not independent transitions, and therefore you cannot simply treat them like we have in the previous models. However, there are two independent transitions: the transition out of \(S\), and the branching to \(A\) or \(I\).

Thus, at each tick \(\Delta t\), susceptible individuals are checked for infection and are moved out of the susceptible compartment with a probability

\[P_\text{out of S} = \lambda_S \left(\frac{A + I}{N}\right)\Delta t.\]

Now, once they are set to transition, they are either sent to \(A\) with a probability \(\gamma\), or otherwise they are sent to \(I\). The asymptomatic and symptomatic individuals are finally transferred to the Removed compartment with a probabilities \(\lambda_A\Delta t\) and \(\lambda_I\Delta t\) respectively.

Once again, we can see the differential equation solutions as the average of the stochastic ones, as demonstrated in the figure below.

The stochastic solutions to the differential equations of the SAIR Model

\(\lambda_S=3/35 \text{ day}^{-1}\) , \(\lambda_A=1/28 \text{ day}^{-1}\), \(\lambda_I=1/28 \text{ day}^{-1}\), \(\gamma=0.6\), \(N=10,000\)

Initial conditions: 1% Asymptomatic, 99% Susceptible

We can now add one last level of complexity to this problem: what if we wanted to model a situation in which asymptomatic individuals are less likely to infect susceptibles (perhaps because they have a lower viral load) than symptomatics. In this case, we would like to include a sort of “relative risk” of infection from an asymptomatic individual that is smaller than the risk of being infected by a symptomatic individual. In order to do this, we can introduce some “contact parameters” that modulate the \(S\to A\) and \(S\to I\) transitions. In this case the differential equations can be written as:

\[\begin{split}\begin{aligned} \dv{S}{t} &= -\frac{\lambda_S}{N} S \left(C_A A + C_I I\right) \\[10pt] \dv{A}{t} &= \gamma \frac{\lambda_S}{N} S\left(C_A A + C_I I\right) - \lambda_A A \\[10pt] \dv{I}{t} &= (1-\gamma) \frac{\lambda_S}{N} S \left(C_A A + C_I I\right) - \lambda_I I \\[10pt] \dv{R}{t} &= \lambda_A A+ \lambda_I I \end{aligned}\end{split}\]

Thus, if \(C_I = 1\) and \(C_A = 0.5\), then a single asymptomatic individual is only half as likely as a symptomatic individual at infecting a susceptible person.

Note

Notice how the quantities that really matter re not \(C_A\) or \(C_I\), but rather \(\lambda_S\, C_A\) and \(\lambda_S\, C_I\). If you were to choose \(C_I = 2\) and \(C_A = 1\), in this case as well asymptomatics will be half as likely like to infect susceptibles, but we have effectively increased the overall value of \(\lambda_S\) because of the factor 2.

Exercise

In this case, would setting \(\lambda_A = \lambda_I\) reduce this to a simple SIR model, as before? Why not?

The SIRS Model

In the SIR model, recovered individuals attain life long immunity. However, this is not the case for many diseases. The acquired immunity can decline over time and as a result the recovered individuals can get reinfected. The SIRS (SusceptibleInfectedRecoveredSusceptible) model allows us to transfer recovered individuals back to the Susceptible compartment. The diagram below shows the movement of the individuals through each compartment in an SIRS model.

The disease progression in the SIRS Model

The infectious rate, \(\lambda_S\), represents the probability of transmitting disease between a susceptible and an infectious individual. \(\lambda_I\) is the recovery rate which can be determined from the average duration of infection.

\(\lambda_R\) is the rate at which the recovered individuals return to the susceptible statue due to loss of immunity.

Ignoring the vital dynamics (births and deaths), in the deterministic form, the SIRS model can be written as the following ordinary differential equations:

\[\begin{split}\begin{aligned} \dv{S}{t} &= -\lambda_S \frac{SI}{N} + \lambda_R R \\[10pt] \dv{I}{t} &= \lambda_S \frac{SI}{N} - \lambda_I I \\[10pt] \dv{R}{t} &= \lambda_I I - \lambda_R R \end{aligned}\end{split}\]

where the total population is,

\[N = S + I + R\]

The main difference between this model and the SIR model is now that because of the possibility of reinfections, there also exists the possibility of multiple _waves_ of infection. In the example below, we can see the emergence of a second wave (easily visible by seeing an increase in \(R(t)\) from days 150-200):

The stochastic solutions to the differential equations of the SIRS Model, demonstrating oscillations

Parameters: \(\lambda_S=2/5 \text{ (day)}^{-1}\), \(\lambda_I=1/5 \text{ (day)}^{-1}\), \(\lambda_R=1/100 \text{ (day)}^{-1}\), \(N=10,000\).

Initial conditions: 1% Infected, 99% Susceptible

On choosing the right parameters, an endemic equilibrium can be reached, meaning that the disease never truly dies out: some small fraction of the population is always infected, as shown below.

The stochastic solutions to the differential equations of the SIRS Model, demonstrating a steady state

Parameters: \(\lambda_S=1/5 \text{ day}^{-1}\), \(\lambda_I=1/20 \text{ day}^{-1}\), \(\lambda_R=1/100 \text{ day}^{-1}\), \(N=10,000\).

Initial conditions: 1% Infected, 99% Susceptible

In the algorithm, during each time step \(\Delta t\), the individuals are transferred from Susceptible to the Infected and from Infected to the Recovered compartments with the same probability as in an SIR model.

\[\begin{split}\begin{aligned} \ P_\text{SI} &= \lambda_S \frac{I}{N} \Delta t\\[10pt] \ P_\text{IR} &= \lambda_I \Delta t \end{aligned}\end{split}\]

The recovered individuals upon loss of immunity are transferred back to the Susceptible compartment with probability,

\[P_\text{RS} = \lambda_R \Delta t\]

Exercises

Basic Exercises

  1. Generalise the simple SIR model to an SEIR model, where the susceptibles move into an “exposed” compartment for some time before they start being infectious. In this compartment, they cannot infect others, nor can they get infected.

  2. Generalise your SIR model to an SIRD model, where some fraction of the infected individuals (say, 1%) transit to the “Dead” compartment instead of recovering.

  3. Generalise the simple SIR model to an SAIR model: infected individuals can either be Asymptomatic (\(A\)) or Symptomatic (\(I\)), and they can both recover from this state. You will need to define new parameters that determine the relative fraction of asymptomatics, for example.

  4. For the simple SIR model (or any other model that you prefer), compute out the “residence time” distribution for the susceptible compartment. In other words, find all the times at which individual agents transition from \(S\to I\), and plot a histogram.

Intermediate Exercises

  1. How does our simple SIR model scale up with the population \(N\)? You will need to study not only the epidemic curves, but also their spread. Run the SIR model for around 10 or 20 times for a given population size (say, 10,000 individuals), and record the output. Repeat the process for a population of 20,000, and so on, for as high as you can go. Next, find some way to quantify the spread of the epidemic curves (for example, the standard deviation in the number of recovered at the end of the epidemic). Plot a graph that shows how this “spread” varies with the initial population \(N\).

  2. Generalise the simple SIR model to include reinfections (the SIRS model): the recovered people do not stay recovered but – at some rate \(\zeta\) – transit back into the susceptible compartment. Try to choose the parameters well such that an endemic equilibrium is reached, meaning that the disease never truly dies out, some small fraction of the population is always infected.

  3. Introduce a lockdown to a simple SIR model. Say that only some fraction of the population (say, only the essential workers) will “violate” the lockdown and continue to go to work. Observe what happens to the total number of recovered at the end of the epidemic as you increase the duration of the lockdown.

Advanced Exercises

  1. Consider a simple model with multiple types of infected individuals (say, the SAIR model I described above). Introduce the possibility that one infected group is much more infectious than another infected group. In other words, say that Symptomatic individuals (“I”) are 2 times more infectious than Asymptomatic individuals.

  2. For the simple SIR model, in one of the earlier exercises we showed that the time that an individual spent in the infected compartment was an exponentially distributed Random Number with mean \(1/\lambda_I\). Therefore, there should be another way to decide how long infected individuals stay in a compartment, instead of generating a probability of transitioning at every step \(\Delta t\). When an individual transits from \(S\to I\), choose a fixed amount of time randomly drawn from an exponential distribution of mean \(\tau_I = 1/\lambda_I\). After this time, the individuals will exit from \(I\to R\).

  3. Compute the effective reproductive ratio (\(R\)) as a function of time. This number is defined as the average number of individuals a single infected individual is responsible for having infected. In order to do this, whenever an \(S\to R\) transition occurs, choose one of the infected individuals who could be responsible for this an increment some parameter representing the the “number of people they have infected”. At every time step, average this number out over the entire population, and finally plot a graph of this as a function of time.

Getting Started

The BharatSim framework is written in Scala 2, so once the source code is obtained, a development environment needs to be set up. In this section, we will describe how to this can be done.

Setup Requirements

BharatSim requires:

  • A JDK (Java Development Kit)

  • A Scala Compiler and SBT (Scala Build Tool)

Installing JDK

BharatSim requires a version of Java that is less than Java 17. We recommend using either Java 8 or Java 11, either of which can be obtained from Oracle or OpenJDK.

Installing Scala with Coursier

Scala can be easily set up using coursier, the Scala application manager. The installation instructions for coursier can be found on the coursier site.

Once the cs binary is downloaded, run:

$ ./cs setup

To check whether the set up was successful, run:

$ scala -version

If the version number is not returned, reboot the system and try running the above line of code again.

Setting up an IDE and First Run

The recommended IDE for BharatSim is IntelliJ Idea. It provides many features helpful for newcomers to a language. It also handles large projects well.

However, the additional features of IntelliJ Idea make it slightly resource intensive. A faster alternative is Microsoft’s Visual Studio Code.

Here we will describe the setup process using both the IntelliJ Idea and VSCode IDEs. However, please note that the BharatSim code can just as well be edited using a simple text editor, and run with the Scala command line sbt tool, described further below. However, this method is not recommended for beginners.

IntelliJ Idea

Note

On Linux, a .tar.gz file is downloaded.

  • Unzip the archive: tar -xvf archivename.tar.gz.

  • Go the unarchived folder. Run the bin/idea.sh file to open IntelliJ Idea.

  • Go to File –> Settings [Or press CTRL+ALT+S] to open the Settings Menu. Go to Plugins and search and install the Scala Plugin.

    _images/idea-scala-plugin.png

    The IDE is now setup.

    Let’s open the Project. The following steps are assuming the source code directory is named BharatSim.

  • Go to File –> Open, and select the build.sbt file in the BharatSim source code directory. Now, select Open as Project. IntelliJ Idea will load the project.

  • Now, we will run the example The SIR Model that is included with the BharatSim source code.

    • In IntelliJ Idea, on the left is the project directory structure and tree. Navigate to BharatSim/src/main/scala/com/bharatsim/examples/epidemiology/sir and open the Main.scala file.

    • There will again be a Green arrow/triangle besides the line containing object Main. Click on it, and Run 'Main'.

      Error

      If it gives an error like Ingestion Failed : java.nio.file.NoSuchFileException: citizen10k.csv, make sure there is a file named citizen10k.csv inside the BharatSim folder. If it is not present, it might have been mistakenly deleted or misplaced. Get the source code again in that case.

    • Wait till the program finishes running. At the end, it should look like this:

      _images/sir-run.png
    • The output CSV file is present at BharatSim/src/main/resources/output_unixtimestamp.csv. This contains the output as per the specification in the program. This can be used to further analyze the results of the SIR Model run.

If you reached till here, Congratulations! IntelliJ Idea is setup and working correctly.

Tip

If the green “Run” arrows do not appear, or some other issues occur and the program does not start to run due to failed dependencies or Scala versions, then there is a simple trick to try.

IntelliJ Idea stores its own configuration for the project inside a .idea directory in the project folder. Delete this folder, and start from scratch, by opening the build.sbt file again and then importing the project.

Visual Studio Code

  • Download the Visual Studio Code. Open VSCode.

  • Go to View –> Extensions and search “Scala”. Install the Scala Syntax (official) and Scala (Metals) extensions.

    _images/vscode-extensions.png

    Let’s open the Project. The following steps are assuming the source code directory is named BharatSim.

  • Go to File –> Open Folder, and select the BharatSim folder. When prompted by VSCode, click on Import Build. This uses an open source tool called sbt to compile and test Scala projects.

    • If you miss it somehow, go to View –> Command Palette [or press CTRL+SHIFT+P] and search for “Import build”. Click on “Metals: Import build” and sit back for a while as VSCode goes through the project structure and builds the project. If you are unable to find such an option, make sure you installed the Metals extension. Restart VSCode if needed.

      Error

      If there is an error notification during the import build process, click on the “more information” option. A new tab will open called Metal Doctor and it will display the source of the error. If the error is in Debugging, then the warning can be ignored and set up process can be carried on.

  • Now, we run the SIR Model. SIR is a simple compartmental model to analyze epidemics, where a person can be either Susceptible (S), Infected (I) or Recovered (R). We will see SIR Model in detail in the Epidemiology section.

    • In VSCode, on the left is the project directory structure and tree. Navigate to BharatSim/src/main/scala/com/bharatsim/examples/epidemiology/sir and open the Main.scala file.

    • There will again be a run | debug above the line containing object Main. Click on run.

      Error

      If it gives an error like Ingestion Failed : java.nio.file.NoSuchFileException: citizen10k.csv, make sure there is a file named citizen10k.csv inside the BharatSim folder. If it is not present, it might have been mistakenly deleted or misplaced. Get the source code again in that case. If the problem persists, then copy the citizen.csv and place it in the main folder BharatSim.

    • Wait till the program finishes running. At the end, it should look like this:

      _images/vscode-sir-run.png
    • The output CSV file is present at BharatSim/src/main/resources/output_unixtimestamp.csv. This contains the output as per the specification in the program. This can be used to further analyze the results of the SIR Model run.

If you reached till here, Congratulations! VSCode is setup and working correctly.

Running Scala on Command Line

Let’s assume the source code directory is named BharatSim. Navigate to the directory in terminal. The sbt tool is often utilized to build a project, which is nothing but compiling, running and testing the project. It also offers the capability of executing each of these processes individually.

  • Compile the project:

    $ sbt compile
    
  • Now, we run the SIR Model. SIR is a simple compartmental model to analyze epidemics, where a person can be either Susceptible (S), Infected (I) or Recovered (R). We will see SIR Model in detail in the Epidemiology section.

    • Do sbt run and wait for a still of main classes to appear on the screen. Select the class number associated to com.bharatsim.examples.epidemiology.sir.Main. It may appear as if the class number is not being typed, but it is! Just input the number and press ENTER . It should start running the simulation.

      Error

      If it gives an error like Ingestion Failed : java.nio.file.NoSuchFileException: citizen10k.csv, make sure there is a file named citizen10k.csv inside the BharatSim folder. If it is not present, it might have been mistakenly deleted or misplaced. Get the source code again in that case. If the problem persists, then copy the citizen.csv and place it in the main folder BharatSim.

      It should look like this:

      _images/cli-sir-run.png
  • The output CSV file is present at BharatSim/src/main/resources/output_unixtimestamp.csv. This contains the output as per the specification in the program. This can be used to further analyze the results of the SIR Model run.

This is how Scala programs can be run through the command line.

Tip

Another way to operate Scala through the command line is to simply type sbt and run the sbt console. The other commands can now be run in succession simply as compile, run and more.

_images/sbt-console.png

Framework Basics

The following sections describe the basics of the BharatSim simulation framework: the different components of a simulation, how to read and write CSV files, how to introduce social interventions, and how to optimise your code for faster performance.

Components of the Simulation Engine

The following sections will explain the basic components and terminology of the BharatSim framework.

Graphs, Nodes, and Relations

_images/FSM_network.png

A sample Graph consisting of different nodes: people, houses, and workplaces. People are connected to their homes and workplaces through relations, and forming a contact network.

The basic data structure in which BharatSim stores its data is a Graph. This Graph is composed of multiple nodes, each specified by a unique 64-bit id. The nodes on our graph are modelled as framework-defined extensions of the Node class. Each of these instances of the Node class can be connected to every other through a Relation.

Error

Only one a single Relation may exist between two instances of a Node class. If multiple relations are defined between two Nodes, then the simulation will not know which Relation to pick, and and error will be thrown.

Agents and network locations are both modelled as extensions of the Node class, while relations are defined between them using the framework-defined addRelation() function to create an underlying network structure. The simulation engine further defines multiple functions that allow the modeler to create, access, update, and delete the nodes and their relations. We will describe them in more detail below.

_images/FSM_class.png

A modeler can define different extensions of the Node class to represent, for example, a Person, a Home, or a Work location.

Agents and Behaviours

In agent-based modelling, a system generally consists of a group of automata that make decisions at every time-step, based on data from each other and the environment. These automata are called “agents”. In BharatSim, agents can be modelled using the framework-defined Agent class, which is an extension of the Node class. To allow for heterogeneity present in real-world individuals, different instances of the Agent class can possess different user-defined attributes, like their age, occupation, vaccination status, and so on. Multiple classes of agents can be defined in a single simulation. Every class of agent will have to be registered using the framework-defined registerAgent() function.

At each time step, an Agent is allowed to execute an action, known as a behaviour. This is implemented in BharatSim through a framework-defined addBehaviour() function that every Agent posesses, which can be used to define a custom behaviour that is executed at every time step.

These behaviours allow the agents in our simulation to mimic the actions of real individuals in a population. For example, in the case of disease-modelling, one might use a behaviour to decide if an unvaccinated agent will get vaccinated on a specific day, based on the result of a daily coin-toss.

Depending on the level of heterogeneity introduced in the population by the modeller, these behaviours can be modelled as close to real-world actions as possible.

_images/FSM_Person.png

Representation of the Agent class: agents can be made to execute actions at every time-step through the addBehaviour function, but can also have custom functions that can be called. The agent’s movement between different network locations is governed by a “schedule” (described below).

Networks

The Network class is another framework-defined extension of the Node class which can be used to model physical locations or contact-networks in a simulation.

In addition to the standard functions that the Node class provides, the Network class has a getContactProbability function which allows the programmer to model differential disease transmission based on network locations. For example, a crowded public-transport location might lead to a much higher probability of transmission of an infectious disease, when compared to an open office with very few employees.

The Network class can be extended by the modeller to describe different network locations, such as homes, workplaces, and schools, for example. Each network location is linked to a class of agent using a different relation, specified by the addRelation function. In this case, the nodes would be the Agent and the specific extension of the Network class.

Note

The addRelation function establishes a relation between any two nodes, not just Agent and Network classes. It is possible to have two agents that have a relation between them, and similarly it is possible to have two network locations that have a relation between them. For an example of the latter, see Houses on a lattice.

_images/FSM_relations.png

Illustration of bidirectional relationships between classes of nodes.

To illustrate the point, consider a simple model in which we define three types of network locations: a Home, a Work location, and a School. Every agent in the population is assigned one of these, with multiple agents being assigned the same home, workplace, and school, based on data from the synthetic population. Each agent is connected to each location using user-defined relations. For example, we could say that a Person IS_EMPLOYED_BY a specific Work location, and that the Work location EMPLOYS the Person. Thus, the relations IS_EMPLOYED_BY and EMPLOYS connect the Person and Work classes.

Every Agent has a Home with a unique id, and therefore a contact network associated with their family – i.e., the other agents who have been assigned the same Home. Similarly, this agent and all other agents who share the same workplace id are assumed to work together, forming a professional network. Agents can be made to spend different amount of times with these different network locations based on their schedules, which can lead to complex social dynamics.

Schedules

In order to account for the movement of individuals between different network locations, the BharatSim framework implements a “schedule” which specifies for how long an Agent is linked to a specific Network location. These schedules allow the modeller to decide what fraction of a unit of time (say, a day of 24 hours) the agent spends in each location. Each agent is assumed to follow this schedule, and this allows them to move between network locations over time, governed by their Schedule.

Different agents can be given different schedules, based on their attributes or other factors. For example, if we define all individuals under the age of 18 to be schoolchildren and all those above the age of 18 to be working adults, we could assign all children to schools and all adults to workplaces, and give them different schedules based on this distinction. Taking another example, a modeller could assign agents an attribute is_employed which is set to true if the agent is employed. We could then define different schedules, an “employed” schedule and an “unemployed” schedule, requiring agents to follow the appropriate schedule based on their employment status. The condition for assigning an Agent a specific schedule can thus be made as general or specific as required.

_images/FSM_customSch.png

Different types of Persons can have their own schedule based on age, jobs, and socio-economic status.

Furthermore, the same agent can be assigned multiple schedules. Which schedule the agent follows at any given simulation tick is dictated by a “priority” parameter that is set when the schedule is defined. With all else being the same, the schedule with the higher priority is given precedence.

Finite State Machine

A Finite State Machine is a class of algorithms where an abstract machine can be in exactly one of the finite state at any given time. A “state” is defined as the explicit trait of the system and this can be changed after satisfying a said boolean condition. This change from one state to another state is called a transition, and the criteria for a transition between two different pair of states may vary based on modeller choices. This is best illustrated through an example of traffic lights.

Traffic Lights as a Finite State Machine

Current State

Next State

Condition

Green

Yellow

120 seconds

Yellow

Red

20 seconds

Red

Green

120 seconds

The above table lists the “state” the system can be in and the possible transition condition that needs to be satisfied. Suppose the system just entered the Green state, then this implies that are after spending 120 seconds being Green, the system will transition to the Yellow state.

BharatSim possesses a framework-defined Finite State Machine which allows users to define a State as an extension of the Scala programming language’s trait datatype. Every distinct state in the model can be created as user-defined extensions of the State trait. These user-defined extensions can further allow for transitions between different states using the framework-defined addTransition function that every State` possesses.

Every type of State should be registered in the Simulation using the framework-defined registerState function.

Stateful Agent

In order to make use of the Finite State Machine, BharatSim defines an extension of the Agent class called the StatefulAgent class. Such agents possess an activeState which links them to one – and only one – instance of a State trait. Additionally, a StatefulAgent also possesses two important functions:

  • setInitialState: which sets the initial state of a StatefulAgent, and

  • fetchActiveState: which returns the current state of the StatefulAgent.

Like agents, a StatefulAgent also needs to be registered in the Simulation using the registerAgent function.

Actions

Additionally, in BharatSim, the State trait also posses certain “actions” that are executed by every StatefulAgent associated with that State, based on certain conditions:

  • An enterAction is an action a StatefulAgent executes the moment they transition into a State for the first time.

  • A perTickAction: is an action executed by a StatefulAgent on every day that they are associated with that State, very much like a behaviour.

Transitions

Additionally, similar to how agents have an addBehaviour function,the State trait allows for an addTransition function that can be used to check at every time-step if each StatefulAgent associated with that State is allowed to transition to another state.

Inputs and Outputs in BharatSim

Reading inputs from a synthetic population

BharatSim uses a CSV file as an input. It is equipped to ingest data from a file, by reading it and converting the data into the Graph used by the framework.

In order to ingest data, we need to make use of three functions:

  • ingestData: which is a method of the Simulation class

  • ingestCSVData: a method of the ContextBuilder object

  • A user-defined “mapper” function which tells the program what data to extract and what to do with it

Using the framework-defined functions

First, we need to import the necessary packages:

import com.bharatsim.engine.ContextBuilder._
import com.bharatsim.engine.execution.Simulation
import com.bharatsim.engine.graph.ingestion.{GraphData, Relation}

Hint

You’ll see why we import GraphData and Relation in the next section!

The next step is to create an instance of the simulation class in the main function,

val simulation = Simulation()

We then ingest the data in the following way:

simulation.ingestData(implicit context => {
  ingestCSVData("input.csv", myCsvDataExtractor)
  logger.debug("Ingestion done")
})

where myCsvDataExtractor is the user-defined “mapper” function.

Note

The above block of code essentially causes the data from the CSV file to be read one line at a time

Creating the user-defined mapping function

The user-defined function (myCsvDataExtractor, in our case) will depend on the data we want to extract. As an example, let’s consider that we have data on a number of cats, each with their own ID, name, city of residence, an integer ID for the city, and a particular colour. Our CSV file would look like

input.csv

ID

Name

City

CityID

Colour

0

Mittens

Sydney

2000

White

1

Tabby

London

1050

Brown

2

Garfield

LasagnaLand

7

Orange

3

Elizabeth

Mishelam

102

Black

4

Coppe

Crossbell

100

Black

5

Marie

Crossbell

100

Orange

6

Antoine

Zeiss

62

Brown

Let’s assume we’ve already defined the following:

  • a case class Cat with three attributes, id, name and colour

  • another case class City with two attributes, id and cityname

Our function myCsvDataExtractor should do the following

  • accept a map of keys (the CSV headers) and values (the CSV element in the row corresponding to the header)

  • accept the context as an implicit parameter

  • return a GraphData object (which you can read about here)

Note

The map is already provided by the ingestCsvData function

private def myCsvDataExtractor(map: Map[String, String](implicit context: Context): GraphData = {}

The first thing we need to do in the function is store the CSV data to appropriate variables.

val catName = map("Name").toString
val catID = map("ID").toLong
val catCity = map("City").toString
val catCityID = map("CityID").toLong
val catColour = map("Colour").toString

Note

The key of the map is the header from the CSV file.

We then use a Constructor to create an instance of the Cat class, for the cat pertaining to a particular row in the CSV. We then do the same for the City class.

val singleCat: Cat = Cat(
  catID,
  catName,
  catColour
)

val singleCity: City = City(
  catCityId,
  catCity
)

Next, we establish relations that will link nodes on the graph. We make a livesIn relation between the cat and the city, and a contains relation between the city and the cat. To do this, we specify the classes the relation is formed between, and then the unique IDs of the nodes with the relation in between them.

val livesIn = Relation[Cat, City](catID, "LIVES_IN", catCityID)
val contains = Relation[City, Cat](catCityID, "CONTAINS", catID)

We then create an instance of the GraphData class, and add the nodes and relations to it

val graphData = GraphData()
graphData.addNode(catID, singleCat)
graphData.addNode(catCityID, singleCity)
graphData.addRelations(staysAt, contains)

Note

The first parameter of graphData.addNode is the unique key of the node.

Finally, we need our function to return the graphData object we’ve made:

graphData

Hint

In Scala, the last line of a function is treated as a return, and so this is valid syntax.

Putting it all together, our user-defined myCsvDataExtractor function is

private def myCsvDataExtractor(map: Map[String, String])(implicit context: Context): GraphData = {

  val catName = map("Name").toString
  val catID = map("ID").toLong
  val catCity = map("City").toString
  val catCityID = map("CityID").toLong
  val catColour = map("Colour").toString

  val singleCat: Cat = Cat(
    catID,
    catName,
    catColour
  )

  val singleCity: City = City(
    catCityId,
    catCity
  )

  val livesIn = Relation[Cat, City](catID, "LIVES_IN", catCityID)
  val contains = Relation[City, Cat](catCityID, "CONTAINS", catID)

  val graphData = GraphData()
  graphData.addNode(catID, singleCat)
  graphData.addNode(catCityID, singleCity)
  graphData.addRelations(staysAt, contains)

  graphData
}

Note

You may have noticed that in the CSV file, two cats (namely, Coppe and Marie) both live in the same city (Crossbell). That does not, however, lead to two nodes being created for the same city. A node is defined by it’s unique key and it’s instance. In this example, the unique key is the city ID (which is the same for both cats - 100) and the instance is the corresponding object singleCity, which is again identical for both the cats (the attributes are 100 and "Crossbell", respectively). As such, the same node is used, and the city doesn’t duplicate in the graph.

Writing outputs to a CSV file

A convenient way to store the output is by using a CSV file. Scala is capable of writing to files, but BharatSim simplifies the process when it comes to CSV outputs.

Note

In case the quantities you’d like to output are fairly simple, you could use Scala’s println function to directly output what you need.

Saving your output to a CSV file

BharatSim relies on a trait called SimulationListener to help output data .

SimulationListener contains 4 methods, each of which allow us to perform a task in one of the following situations:

  • At the start of the simulation

  • At the start of every time step

  • At the end of every time step

  • At the end of the simulation

The BharatSim engine also contains a class called CsvOutputGenerator, an extension of SimulationListener which has two attributes:

  • path, the desired path for the output file to be stored

  • csvSpecs, a user-defined class that outputs the headers and the rows required. Note that this user-defined class should extend the CSVSpecs trait and override the getHeaders and getRows methods.

This class writes the headers at the start of the simulation, writes the rows at the start of every time step, and closes the writer at the end of the simulation.

Output at a single instant of time

We can define a class as follows:

import com.bharatsim.engine.Context
import com.bharatsim.engine.listeners.CSVSpecs

class MyOutputSpec(context: Context) extends CSVSpecs {
  override def getHeaders: List[String] =
    List(
      "Header1",
      "Header2",
      "Header3"
    )
  override def getRows(): List[List[Any]] = {
    val elementInRow: String = "row" + context.getCurrentStep.toString
    val row = List(
      elementInRow,
      elementInRow,
      elementInRow
    )
    List(row)
  }
}

Now, we need to create an instance of the CsvOutputGenerator class that uses MyOutputSpec, and call the required methods. First, we need to import CsvOutputGenerator into our main class:

import com.bharatsim.engine.listeners.CsvOutputGenerator

Next, we add the following code snippet inside simulation.defineSimulation in the main function:

var outputGenerator = new CsvOutputGenerator("src/main/resources/output.csv", new MyOutputSpec(context))
outputGenerator.onSimulationStart(context)
outputGenerator.onStepStart(context)
outputGenerator.onSimulationEnd(context)

Note

Calling the onStepEnd method of the class isn’t necessary, as the CsvOutputGenerator class currently does nothing when it’s called.

The output is

output.csv

Header1

Header2

Header3

row0

row0

row0

Hint

In case you want your outputs generated after the simulation is completed, you can place the above 4 lines of code inside simulation.onCompleteSimulation.

You can see a more in-depth example of this in Saving location-level information from the simulation.

Output at every time step

If we’d like to investigate the dynamics of the simulation as it evolves with time, we essentially need to call the three methods described above every time step. BharatSim simplifies things with SimulationListenerRegistry, which allows us to register the output generator in the simulation (similar to how we registered agents), so that it writes data to the CSV file at every time step.

First, we must import CsvOutputGenerator and SimulationListenerRegistry

import com.bharatsim.engine.listeners.{CsvOutputGenerator, SimulationListenerRegistry}

Next, we register it using the register method of SimulationListenerRegistry. Note that the following code snippet must go inside simulation.defineSimulation in the main function.

SimulationListenerRegistry.register(
  new CsvOutputGenerator("src/main/resources/output.csv", new myOutputSpec(context))
  )

where myCsvSpecs is the user-defined class which requires the context as an attribute.

Now, the output is

output.csv

Header1

Header2

Header3

row1

row1

row1

row2

row2

row2

row3

row3

row3

row4

row4

row4

row5

row5

row5

and so on, until the tick at which the simulation ends.

Hint

Running the above block of code once will cause a file called output to be created at src/main/resources/. However, running it again will rewrite the contents of the file with the new output. You can get around this by adding the current time to the output as a string. For example,

val currentTime = new Date().getTime

SimulationListenerRegistry.register(
    new CsvOutputGenerator("src/main/resources/output_" + currentTime + ".csv", new SIROutputSpec(context))
  )

Note that Date().getTime returns the time as a UNIX timestamp, and so your output will contain a long integer after the underscore.

For a more detailed example of how to output data to a CSV file, please refer to the Writing your first program section.

Interventions in BharatSim

Interventions are punctual events that occur when some external condition is satisfied. These could, for example, be policy-level decisions made by governments, or environmental factors like the introduction of new viral strains, for example. Each intervention is identified uniquely by the the Name of the intervention. Additionally, each intervention requies an activation condition, and a deactivation condition. These conditions are boolean decisions.

Additionally, the user can further define an optional activationAction and a perTickAction.

  • activationAction: an action that is invoked whenever the intervention is activated.

  • perTickAction: an action that is invoked for every tick for which the intervention is active.

There are four different intervention classes available in BharatSim

Intervention

This is a generic intervention that can be invoked whenever the activation condition is met based on the parameters passed. This can be used to implement lockdowns during an epidemic every time the infected fraction of agent reaches a certain threshold, for example.

Five parameters can be passed Intervention object:

  • interventionName: This is the unique name of the intervention.

  • activationCondition: This function tells whether this intervention should be activated or not. The function should return a Boolean value.

  • deActivationCondition: This function tells whether this intervention should be deactivated. The function should return a Boolean value.

  • firstTimeExecution: This is an optional function that will be executed at the start of the intervention

  • whenActiveActionFunc: This is an optional function executed per tick when intervention is active.

Example using the Intervention object

Here we will use the Intervention class to implement a two-week lockdown every time the number of infected agents are greater than or equal to 2000. First, we will define the intervention named addLockdown and initialise the variables.

interventionActivatedAt stores the information about when the intervention was activated. Here it is initialised to be zero and will be updated once the intervention is activated.

activationCondition has a boolean value that defines when the intervention has to be activated. Here the intervention gets activated once the number of infected agents is greater than or equal to 2000.

firstTimeExecution specifies what should be done once the intervention is activated. This is only executed once when the intervention is activated

deActivationCondition specifies the condition when the intervention should be stoped. Here it is stopped after 14 days, i.e. 2*14 ticks after the intervention activation.

  private def addLockdown(implicit context: Context): Unit = {

  var interventionActivatedAt = 0
  val interventionName = "lockdown"
  val activationCondition = (context: Context) => getInfectedCount(context) >= 2000
  val firstTimeExecution = (context: Context) => interventionActivatedAt = context.getCurrentStep
  val deActivationCondition = (context: Context) => {
  context.getCurrentStep >= interventionActivatedAt + 2*14
  }
}

Now we have to create an instance of the Intervention object and pass the values defined earlier it:

val intervention =
Intervention(interventionName, activationCondition, deActivationCondition, firstTimeExecution)

Next, we define a new schedule that has to come into effect once the intervention is activated. Here we define the lockdownSchedule such that all agent stays in the house throughout the day. The code below makes the agent stay home for the entire day as long as the intervention remains activated.

val lockdownSchedule = (myDay, myTick).add[House](0, 1)

Now we have to register the register the intervention and the schedules using registerIntervention and registerSchedules respectively. We also have to pass the Agent and Context to registerSchedules.

registerIntervention(intervention)
registerSchedules(
  (
    lockdownSchedule,
    (agent: Agent, context: Context)
  )
)

The complete definition of our user-defined addLockdown is given below.

private def addLockdown(implicit context: Context): Unit = {

var interventionActivatedAt = 0
val interventionName = "lockdown"
val activationCondition = (context: Context) => getInfectedCount(context) >= 2000
val firstTimeExecution = (context: Context) => interventionActivatedAt = context.getCurrentStep
val deActivationCondition = (context: Context) => {
context.getCurrentStep >= interventionActivatedAt + 2*14
}
val intervention =
Intervention(interventionName, activationCondition, deActivationCondition, firstTimeExecution)

val lockdownSchedule = (myDay, myTick).add[House](0, 1)

registerIntervention(intervention)
registerSchedules(
(
    lockdownSchedule,
    (agent: Agent, context: Context)
)
)
}

Hint

addLockdown should be called in the definition of the simulation.

simulation.defineSimulation(implicit context => {
addLockdown
}

IntervalBasedIntervention

This can be used to invoke an intervention that starts and end at specific ticks. This can be used for giving relaxations in the epidemic regulations for a specified period, for example, during a festival.

Five paramertes can be passed to the IntervalBasedIntervention object:

  • interventionName: This is the unique intervention name.

  • startTick: This integer specifies the start tick for intervention (inclusive); it should not be greater than endTick.

  • endTick: This is an integer that specifies the end tick for the intervention (It is exclusive, and intervention will not be active at “endTick”.)

  • firstTimeActionFunc: This is an optional function which gets executed when simulation starts.

  • whenActiveActionFunc: This is an optional function executed per tick when the simulation is active.

Example using the IntervalBasedIntervention object

We will try to implement during the 50th to the 100th tick, i.e. for 25 days.

First we will define the variables interventionName. In the function call of IntervalBasedIntervention(), we will pass the interventionName, startTick and endTick.

private def addLockdown(implicit context: Context): Unit = {

val interventionName = "lockdown"
val intervention = IntervalBasedIntervention(interventionName, 50, 100)
}

Now we define the lockdownSchedule to force all agents to stay home for the entire day throughout all ticks when the intervention is active.

val lockdownSchedule = (Day, Hour).add[House](0, 1)

Now we will register both the intervention as well as the schedule.

registerIntervention(intervention)
  registerSchedules(
  (

      lockdownSchedule,
      (agent: Agent, context: Context)

  )
  )

The entire definition of addLockdown intervention using the IntervalBasedIntervention is given below:

private def addLockdown(implicit context: Context): Unit = {

val interventionName = "lockdown"
val intervention = IntervalBasedIntervention(interventionName, 50, 100)

val lockdownSchedule = (Day, Hour).add[House](0, 1)

registerIntervention(intervention)
registerSchedules(
(

    lockdownSchedule,
    (agent: Agent, context: Context)

)
)
}

Hint

addLockdown should be included in the definition of the simulation.

simulation.defineSimulation(implicit context => {
addLockdown
}

OffsetBasedIntervention

This can be used to invoke interventions that end after ‘n’ ticks. It gets invoked when the shouldActivateWhen function is true. This can be used to implement a lockdown when the number of infected agents reaches a particular threshold and stays active till n ticks.

Five parameters can be passed to the OffsetBasedIntervention object:

  • interventionName: This is the unique name of the intervention.

  • shouldActivateWhen: This function decides when should the intervention be activated.

  • endAfterNTicks: This is the offset ‘n’; simulation will end after n ticks from the star tick.

  • firstTimeActionFunc:This is an optional function which gets executed when simulation starts.

  • whenActiveActionFunc: This is an optional function executed per tick when the simulation is active.

Example suing the OffsetBasedIntervention object

We will implement a lockdown when the number of infected individuals is greater than or equal to 2000 and stays active for 28 ticks (14 days) from the start of the lockdown. First, we will define the intervention named addLockdown and initialise the variables.

activationCondition has a boolean value that defines when the intervention has to be activated. Here the intervention gets activated once the number of infected agents interventionName contains the name of the intervention

  private def addLockdown(implicit context: Context): Unit = {

  val interventionName = "lockdown"
  val activationCondition = (context: Context) => getInfectedCount(context) >= 2000
  }
}

Now we have to create an instance of the Intervention object. Here we define intervention and pass the interventionName, activationCondition and the number of ticks after which the intervention has to stop to the Intervention object.

val intervention =
Intervention(interventionName, activationCondition,28)

Now we define a new schedule that has to come into effect once the intervention is activated. Here we define the lockdownSchedule such that all agent stays in the house throughout the day. Here the 0,1 passed to the add[House] makes the agent stay in its house from tick 0 to the end of tick 1 in a day. Here one day is defined as having two ticks, i.e. 0 and 1. So this makes the agent stay home for the entire day as long as the intervention remains activated.

val lockdownSchedule = (myDay, myTick).add[House](0, 1)

Now we have to register the register the intervention and the schedules using registerIntervention and registerSchedules respectively. We also have to pass the Agent and Context to registerSchedules.

registerIntervention(intervention)
registerSchedules(
  (
    lockdownSchedule,
    (agent: Agent, context: Context)
  )
)

The complete definition of addLockdown using OffsetBasedIntervention is given below.

private def addLockdown(implicit context: Context): Unit = {

var interventionActivatedAt = 0
val interventionName = "lockdown"
val activationCondition = (context: Context) => getInfectedCount(context) >= 2000
val firstTimeExecution = (context: Context) => interventionActivatedAt = context.getCurrentStep
val deActivationCondition = (context: Context) => {
context.getCurrentStep >= interventionActivatedAt + 2*14
}
val intervention =
Intervention(interventionName, activationCondition, deActivationCondition, firstTimeExecution)

val lockdownSchedule = (myDay, myTick).add[House](0, 1)

registerIntervention(intervention)
registerSchedules(
(
    lockdownSchedule,
    (agent: Agent, context: Context)
)
)
}

Hint

As before, addLockdown should be called in the definition of the simulation.

simulation.defineSimulation(implicit context => {
addLockdown
}

SingleInvocationIntervention

This can be used to create an intervention that will be invoked only once in the simulation.

Five parameters can be passed to the SingleInvocationIntervention object:

  • interventionName: This is the unique name of the intervention.

  • shouldActivateFunc: This function tells whether this intervention should be activated.

  • shouldDeactivateFunc: This function tells whether this intervention should be deactivated.

  • firstTimeActionFunc: This is an optional function executed at the start of the intervention.

  • whenActiveActionFunc: This is an optional function executed per tick when intervention is active.

Example using the SingleInvocationIntervention object

Here we will use the Intervention class to implement a lockdown just once when the number of infected agent crosses 2000.

First, we will define the intervention named addLockdown and initialise the variables.

interventionActivatedAt stores the information about when the intervention was activated. Here it is initialised to be zero and will be updated once the intervention is activated.

activationCondition has a boolean value that defines when the intervention has to be activated. Here the intervention gets activated once the number of infected agents is greater than or equal to 2000.

firstTimeExecution specifies what should be done once the intervention is activated. This is only executed once when the intervention is activated

deActivationCondition specifies the condition when the intervention should be stoped. Here it is stopped after 14 days, i.e. 2*14 ticks after the intervention activation.

  private def addLockdown(implicit context: Context): Unit = {

  var interventionActivatedAt = 0
  val interventionName = "lockdown"
  val activationCondition = (context: Context) => getInfectedCount(context) >= 2000
  val firstTimeExecution = (context: Context) => interventionActivatedAt = context.getCurrentStep
  val deActivationCondition = (context: Context) => {
  context.getCurrentStep >= interventionActivatedAt + 2*14
  }
}

Now we have to create an instance of the Intervention object. Here we define intervention and pass the values defined in earlier` to the Intervention object.

val intervention =
SingleInvocationIntervention(interventionName, activationCondition, deActivationCondition, firstTimeExecution)

Now we define a new schedule that has to come into effect once the intervention is activated. Here we define the lockdownSchedule such that all agent stays in the house throughout the day. Here the 0,1 passed to the add[House] makes the agent stay in its house from tick 0 to the end of tick 1 in a day. Here one day is defined as having two ticks, i.e. 0 and 1. So this makes the agent stay home for the entire day as long as the intervention remains activated.

val lockdownSchedule = (myDay, myTick).add[House](0, 1)

Now we have to register the register the intervention and the schedules using registerIntervention and registerSchedules respectively. We also have to pass the Agent and Context to registerSchedules.

registerIntervention(intervention)
registerSchedules(
  (
    lockdownSchedule,
    (agent: Agent, context: Context)
  )
)

The complete definition of addLockdown using SingleInvocationIntervention is given below.

private def addLockdown(implicit context: Context): Unit = {

var interventionActivatedAt = 0
val interventionName = "lockdown"
val activationCondition = (context: Context) => getInfectedCount(context) >= 2000
val firstTimeExecution = (context: Context) => interventionActivatedAt = context.getCurrentStep
val deActivationCondition = (context: Context) => {
context.getCurrentStep >= interventionActivatedAt + 2*14
}
val intervention =
SingleInvocationIntervention(interventionName, activationCondition, deActivationCondition, firstTimeExecution)

val lockdownSchedule = (myDay, myTick).add[House](0, 1)

registerIntervention(intervention)
registerSchedules(
(
    lockdownSchedule,
    (agent: Agent, context: Context)
)
)
}

Hint

addLockdown should be included in the definition of the simulation.

simulation.defineSimulation(implicit context => {
addLockdown
}

Writing your First Program

This section is a detailed guide for a novice user on how to build a SIR Model in BharatSim from scratch. Going through the instructions given in these sections, you should be able to write your own SIR model with agents, introduce network structures like homes and workplaces and allow agents to move between them, and define different disease compartments representing different disease states and transition agents between them.

Any model built on BharatSim contains various classes which are essentially different extensions of the Node class. To build a SIR model from scratch, one needs to define these classes and the properties and relationships associated with them. In what follows, we will be modelling the individuals of our population as extensions of the Agent class, and the locations that they frequent as extensions of the Network class.

Note

Before we move on, make sure you have gone through the sections on Agents and Behaviours and Networks, since certain ideas in the sections below are explained in greater detail there.

Different agents come in contact with each other based on the underlying contact network structure of the population. For simplicity, we begin our model assuming that every agent is in contact with every other agent, which is equivalent to all of them being in the same “location”. However, we will relax this very quickly to account for a more realistic contact network strucure. Finally, we will describe BharatSim’s Finite State Machine (FSM) and rewrite our model using it, with the help of the framework-defined StatefulAgent class.

Single Location SIR

This section will look at the disease progression in a single location and observe its dynamics.

Creating an Empty Class

Create a new project in the required directory. There is no need to create a new folder, as creating a new project automatically creates a new folder. Name this folder sir and change the language to scala. The build system should be chosen to be sbt.

Navigate to src\main\scala and right click on the folder in IntelliJ. Select a new package and rename it sir. Again right click on sir package, select a Scala class followed by object. Call this object Main.

The empty class should look like this.

_images/New_class.png

Now we can define a main function that has no input and has no output. The syntax and indentation of defining a function is as follows

def main(args: Array[String]): Unit = {
}

The args means that the argument or the input is an array of Strings and the output is of type Unit, which corresponds to void means that there is no output. The code should look like this,

Note

void return means that the function returns nothing at all. Remember nothing is different from 0, or empty list.

_images/EmptyClass.png

Note

Notice how the object Main has changed color from grey to white. This is an IntelliJ feature which lets the user know if the object/class/variable is being used

On running this, the output message should read Process finished with exit code 0

Implementing a single-location SIR

Before we can input a file or simulate a disease, we need to make a few classes which are essential to the workings of the framework. These classes need to be imported to the main class to make the code easier to understand and clutter-free. The framework is extremely inter-connected and defining the same functions over and over again is tedious and computationally heavy.

InfectionStatus class is a scala object class that stores the compartments of the disease, and in our case Susceptible, Infected, and Recovered. This class connects the instance of the compartments to the their string counterparts.

  • BasicDecoder: This functions takes a string value and converts it to either a node or throws an exception. The latter is only the case when the input type is not in form of a string.

  • BasicEncoder: This takes the instance and converts it to a string. In the simple case, there are three possibilities which are Susceptible, Infected and Recovered

  • Extends: This allows the functions of one class to be used in another. In this case, the functions of Enumeration are made available in the class InfectedStatus because of extends

The code for each of the above class is provided below.

package sir
import com.bharatsim.engine.basicConversions.StringValue
import com.bharatsim.engine.basicConversions.decoders.BasicDecoder
import com.bharatsim.engine.basicConversions.encoders.BasicEncoder

object InfectionStatus extends Enumeration {
  type InfectionStatus = Value
  val Susceptible, Infected, Removed = Value

  implicit val infectionStatusDecoder: BasicDecoder[InfectionStatus] = {
    case StringValue(v) => withName(v)
    case _ => throw new RuntimeException("Infection status was not stored as a string")
  }

  implicit val infectionStatusEncoder: BasicEncoder[InfectionStatus] = {
    case Susceptible => StringValue("Susceptible")
    case Infected => StringValue("Infected")
    case Removed => StringValue("Removed")
  }
}

Inputting a File

To begin we must import a series of libraries and the function of each libraries will be explained as and when they are required.

import com.bharatsim.engine.Context
import com.bharatsim.engine.ContextBuilder._
import com.bharatsim.engine.execution.Simulation
import com.bharatsim.engine.graph.ingestion.{GraphData, Relation}
import com.typesafe.scalalogging.LazyLogging
import com.bharatsim.engine.utils.Probability.biasedCoinToss
import com.bharatsim.engine.basicConversions.encoders.DefaultEncoders._

There needs to be a modification in the line where we have defined the object. We need to make use of a keywork called extends which allows one class to inherit the properties of another class.

object Main extends LazyLogging

By extending LazyLogging, all the properties of this class are made available in Main. The LazyLogging class allows the user to display or output information. It can be thought of as better version of SystemOut.

Note

When libraries or variables are not being used they appear grey in color, and as soon as they are called, they become colored again

Since LazyLogging is being used, it changes color from grey.

The next step is to define a private value called initialInfectedFraction and set it to 0.01. Private value means that this will only be available in the defining class and not outside. This will be made accessible to the function we are about to define.

In the main function we had earlier defined, we can create an instance of the simulation class.

val simulation = Simulation()

Note

val is an immutable variable and this implies that the value of this can not change.

Then we ingest the csv file in the following manner

simulation.ingestData(implicit context => {
ingestCSVData("input.csv", csvDataExtractor)
logger.debug("Ingestion done")
})

Here csvDataExtractor is a user defined function which we will get to later.

On running the code, an error pops up displaying that csvDataExtractor is not defined.

The csvDataExtractor function is defined in the following manner

private def csvDataExtractor(map: Map[String, String])(implicit context: Context): GraphData = {
}

Once the function is defined and we need it to the following things,

  1. Accept the Context as an input parameter

  2. CSV header and corresponding values

  3. Return the data in the form of GraphData

The first step depends on the CSV file that is being imported since it depends on the headers of the data. In BharatSim, the CSV files usually have the following columns,

val citizenId = map("Agent_ID").toLong
val age = map("Age").toInt
val homeId = map("HHID").toLong

Note

The csvDataExtractor reads the csv file line by line and defines each citizen line by line.

The next step is to determine if the citizen imported is infected or not.

val initialInfectionState = if (biasedCoinToss(initialInfectedFraction)) "Infected" else "Susceptible"

If the biasedCoinToss returns True, then the citizen analyzed is infected from the disease. Using the data obtained from the CSV file and the infection state, we can create an instance of the citizen.

val citizen: Person = Person(
citizenId,
age,
InfectionStatus.withName(initialInfectionState),
0
)

Once this is done, relationships need to be established that will connect the nodes on the graph. The citizen will Stay At the house, and the house will House the citizen. The relationship needs to be established both the ways, as the first relationship links the citizen node to the house node and the second one links the house node to the citizen one.

val home = House(homeId)
val staysAt = Relation[Person, House](citizenId, "STAYS_AT", homeId)
val memberOf = Relation[House, Person](homeId, "HOUSES", citizenId)

Note

A House HOUSES an Agent and an Agent STAYS_AT a House so these two relations are inherently reflections of each other. The first relation is specified in the House class, while the second one is specified in the Person class (Refer to the classes above). The same defination of relationships can be extended to any pair of Agents (Student, Employer) and corresponding locations (School, Office).

Then we create an instance of the GraphData and add the aforementioned nodes and relationships

val graphData = GraphData()
graphData.addNode(citizenId, citizen)
graphData.addNode(homeId, home)
graphData.addRelations(staysAt, memberOf)

Once the nodes and relationships have been established, we can then return the GraphData. Unlike python, no return keywork is actually required. In scala, the last line has to be just value that has to be returned.

graphData

Compiling all the lines together, the csvDataExtractor function and the main function looks like

def main(args: Array[String]): Unit = {

  var beforeCount = 0
  val simulation = Simulation()

  simulation.ingestData(implicit context => {
    ingestCSVData("citizen10k.csv", csvDataExtractor)
    logger.debug("Ingestion done")
  })

private def csvDataExtractor(map: Map[String, String])(implicit context: Context): GraphData = {

  val citizenId = map("Agent_ID").toLong
  val age = map("Age").toInt
  val homeId = map("HHID").toLong

  val initialInfectionState = if (biasedCoinToss(initialInfectedFraction)) "Infected" else "Susceptible"

  val citizen: Person = Person(
    citizenId,
    age,
    InfectionStatus.withName(initialInfectionState),
    0
  )

  val home = House(homeId)
  val staysAt = Relation[Person, House](citizenId, "STAYS_AT", homeId)
  val memberOf = Relation[House, Person](homeId, "HOUSES", citizenId)

  val graphData = GraphData()
  graphData.addNode(citizenId, citizen)
  graphData.addNode(homeId, home)
  graphData.addRelations(staysAt, memberOf)

  graphData
}

Introduction of Disease Dynamics

In the previous section, while we defined all of the disease parameters, it had no effect on the population since we did not allow it to spread or die out. In this section, we allow the disease to propagate through a population and we output the changes in the population, such the number of individuals that have been recovered or number of infected individuals that remain after end time. Since the manner in which the disease interacts with an agent can now change as a function of location and time, we will have to update the Person class to account for this. Additionally, since we would like to record the different number of individuals in each compartment as the output for our simulation, we also need to define a new class that can record these details to an output file.

The Required Classes

A new class called SIROutputSpec needs to be created and the Person class needs to be updated.

As mentioned earlier, this is the updated version of the class we have written earlier. In the previous version, we had only defined the relation and nothing else. The first thing to do is to add a schedule followed by checking the InfectedStatus of the individuals and the people around. The latter is done so we can look at the probability of getting infected and then do a coin toss with this probability to determine if the person in question does get infected.

  • numberOfTicksInADay is used to define how many Ticks a person experiences is a day. Since the duration of the infection (in days) is fixed, the numberOfTicksInADay will dictate the increments in the simulation.

  • incrementInfectionDuration updates the day in the simulation. This is done after all the ticks have been completed in the day, and only after this can we move to the next day.

  • checkForInfection is a function that is used to check whether a susceptible individual gets infected. If the location is not empty, then the number of people present at that location are counted and are infected and this is stored as infectedNeighbourCount. Using these value, an appropriate biased coin toss is done and if it comes True, then the susceptible individual contracts the disease. The InfectionStatus will changed from susceptible to infected

  • checkForRecovery looks at infected individuals and if the last day for infection has been reached, then the InfectionStatus changes from Infected to Recovered.

  • isSusceptible, isInfected, isRecovered changes the infection status to Susceptible, Infected, Recovered respectively.

  • decodeNode take the string and return the corresponding node.

  • We then add behaviour for each of the states.

package sir

import com.bharatsim.engine.Context
import com.bharatsim.engine.basicConversions.decoders.DefaultDecoders._
import com.bharatsim.engine.basicConversions.encoders.DefaultEncoders._
import com.bharatsim.engine.graph.GraphNode
import com.bharatsim.engine.models.{Agent, Node}
import com.bharatsim.engine.utils.Probability.toss
import sir.InfectionStatus._

case class Person(id: Long, age: Int, infectionState: InfectionStatus, infectionDay: Int) extends Agent {
  final val numberOfTicksInADay: Int = 24
  private val incrementInfectionDuration: Context => Unit = (context: Context) => {
    if (isInfected && context.getCurrentStep % numberOfTicksInADay == 0) {
      updateParam("infectionDay", infectionDay + 1)
    }
  }
  private val checkForInfection: Context => Unit = (context: Context) => {
    if (isSusceptible) {
      val infectionRate = Disease.beta

      val schedule = context.fetchScheduleFor(this).get

      val currentStep = context.getCurrentStep
      val placeType: String = schedule.getForStep(currentStep)

      val places = getConnections(getRelation(placeType).get).toList

      if (places.nonEmpty) {
        val place = places.head
        val decodedPlace = decodeNode(placeType, place)

        val infectedNeighbourCount = decodedPlace
          .getConnections(decodedPlace.getRelation[Person]().get)
          .count(x => x.as[Person].isInfected)

        val shouldInfect = toss(infectionRate, infectedNeighbourCount)
        if (shouldInfect) {
          updateParam("infectionState", Infected)
        }
      }
    }
  }

  private val checkForRecovery: Context => Unit = (context: Context) => {
    if (isInfected && infectionDay == Disease.lastDay
    )
      updateParam("infectionState", Removed)
  }

  def isSusceptible: Boolean = infectionState == Susceptible

  def isInfected: Boolean = infectionState == Infected

  def isRecovered: Boolean = infectionState == Removed


  private def decodeNode(classType: String, node: GraphNode): Node = {
    classType match {
      case "House" => node.as[House]
    }
  }

  addBehaviour(incrementInfectionDuration)
  addBehaviour(checkForInfection)
  addBehaviour(checkForRecovery)

  addRelation[House]("STAYS_AT")
}

Writing output to a file

Now we have imported a population and set up basics for the disease. It is time we implement the disease and print the output. First we need to import the following addition files,

import sir.InfectionStatus._
import com.bharatsim.engine.{Context, Day, Hour, ScheduleUnit}
import com.bharatsim.engine.actions.StopSimulation
import com.bharatsim.engine.listeners.{CsvOutputGenerator, SimulationListenerRegistry}
import com.bharatsim.engine.models.Agent
import java.util.Date
import com.bharatsim.engine.basicConversions.decoders.DefaultDecoders._
import com.bharatsim.engine.graph.patternMatcher.MatchCondition._
import com.bharatsim.engine.dsl.SyntaxHelpers._

After we ingest the data in the main function, we need to define the Simulation and the end point of the Simulation. registerAgent[Person] explicitly mentions that the individual of the person class is an agent in the system. Once we define the output location, we can actually run the simulation followed by printing the results, and finally saving the data as a csv file.

def main(args: Array[String]): Unit = {

  var beforeCount = 0
  val simulation = Simulation()

  simulation.ingestData(implicit context => {
    ingestCSVData("citizen10k.csv", csvDataExtractor)
    logger.debug("Ingestion done")
  })

  simulation.defineSimulation(implicit context => {

    createSchedules()

    registerAction(
      StopSimulation,
      (c: Context) => {
        getInfectedCount(c) == 0
      }
    )

    beforeCount = getInfectedCount(context)

    registerAgent[Person]

    val currentTime = new Date().getTime

    SimulationListenerRegistry.register(
      new CsvOutputGenerator("src/main" + currentTime + ".csv", new SIROutputSpec(context))
    )
  })

  simulation.onCompleteSimulation { implicit context =>
    printStats(beforeCount)
    teardown()
  }

  val startTime = System.currentTimeMillis()
  simulation.run()
  val endTime = System.currentTimeMillis()
  logger.info("Total time: {} s", (endTime - startTime) / 1000)
}

In the defineSimulation, we call upon a function called createSchedules. The following piece of code will define this function

private def createSchedules()(implicit context: Context): Unit = {
  val allSchedule = (Day, Hour)
    .add[House](0, 23)

  registerSchedules(
    (allSchedule, (agent: Agent, _: Context) => agent.asInstanceOf[Person].age > 0, 1),
  )
}

Note

add[House](0,23) means that we are creating a 24 hour schedule associated with the location House. In the framework, 0 to 0 is counted as 1 hour.

printStats simply prints the values in the output message window and it finds these values by calling user defined like getSusceptibleCount. These functions look at the node on the graph and then count the people present in the node.

private def printStats(beforeCount: Int)(implicit context: Context): Unit = {
  val afterCountSusceptible = getSusceptibleCount(context)
  val afterCountInfected = getInfectedCount(context)
  val afterCountRecovered = getRemovedCount(context)

  logger.info("Infected before: {}", beforeCount)
  logger.info("Infected after: {}", afterCountInfected)
  logger.info("Recovered: {}", afterCountRecovered)
  logger.info("Susceptible: {}", afterCountSusceptible)
}

private def getSusceptibleCount(context: Context) = {
  context.graphProvider.fetchCount("Person", "infectionState" equ Susceptible)
}

private def getInfectedCount(context: Context): Int = {
  context.graphProvider.fetchCount("Person", ("infectionState" equ Infected))
}

private def getRemovedCount(context: Context) = {
  context.graphProvider.fetchCount("Person", "infectionState" equ Removed)
}

On Compiling everything together, the whole code looks like the following

package sir
import com.bharatsim.engine.Context
import com.bharatsim.engine.ContextBuilder._
import com.bharatsim.engine.execution.Simulation
import com.bharatsim.engine.graph.ingestion.{GraphData, Relation}
import com.typesafe.scalalogging.LazyLogging
import com.bharatsim.engine.utils.Probability.biasedCoinToss
import com.bharatsim.engine.basicConversions.encoders.DefaultEncoders._
import sir.InfectionStatus._
import com.bharatsim.engine.{Context, Day, Hour, ScheduleUnit}
import com.bharatsim.engine.actions.StopSimulation
import com.bharatsim.engine.listeners.{CsvOutputGenerator, SimulationListenerRegistry}
import com.bharatsim.engine.models.Agent
import java.util.Date
import com.bharatsim.engine.basicConversions.decoders.DefaultDecoders._
import com.bharatsim.engine.graph.patternMatcher.MatchCondition._
import com.bharatsim.engine.dsl.SyntaxHelpers._

object Main extends LazyLogging{
  private val initialInfectedFraction = 0.01

  def main(args: Array[String]): Unit = {

    var beforeCount = 0
    val simulation = Simulation()

    simulation.ingestData(implicit context => {
      ingestCSVData("citizen10k.csv", csvDataExtractor)
      logger.debug("Ingestion done")
    })
    simulation.defineSimulation(implicit context => {

      createSchedules()

      registerAction(
        StopSimulation,
        (c: Context) => {
          getInfectedCount(c) == 0
        }
      )

      beforeCount = getInfectedCount(context)

      registerAgent[Person]

      val currentTime = new Date().getTime

      SimulationListenerRegistry.register(
        new CsvOutputGenerator("src/main" + currentTime + ".csv", new SIROutputSpec(context))
      )
    })

    simulation.onCompleteSimulation { implicit context =>
      printStats(beforeCount)
      teardown()
    }

    val startTime = System.currentTimeMillis()
    simulation.run()
    val endTime = System.currentTimeMillis()
    logger.info("Total time: {} s", (endTime - startTime) / 1000)
  }

  private def createSchedules()(implicit context: Context): Unit = {
    val allSchedule = (Day, Hour)
      .add[House](0, 23)

    registerSchedules(
      (allSchedule, (agent: Agent, _: Context) => agent.asInstanceOf[Person].age > 0, 1),
    )
  }
  private def csvDataExtractor(map: Map[String, String])(implicit context: Context): GraphData = {

    val citizenId = map("Agent_ID").toLong
    val age = map("Age").toInt
    val homeId = map("HHID").toLong

    val initialInfectionState = if (biasedCoinToss(initialInfectedFraction)) "Infected" else "Susceptible"

    val citizen: Person = Person(
      citizenId,
      age,
      InfectionStatus.withName(initialInfectionState),
      0
    )

    val home = House(homeId)
    val staysAt = Relation[Person, House](citizenId, "STAYS_AT", homeId)
    val memberOf = Relation[House, Person](homeId, "HOUSES", citizenId)

    val graphData = GraphData()
    graphData.addNode(citizenId, citizen)
    graphData.addNode(homeId, home)
    graphData.addRelations(staysAt, memberOf)

    graphData
  }

  private def printStats(beforeCount: Int)(implicit context: Context): Unit = {
    val afterCountSusceptible = getSusceptibleCount(context)
    val afterCountInfected = getInfectedCount(context)
    val afterCountRecovered = getRemovedCount(context)

    logger.info("Infected before: {}", beforeCount)
    logger.info("Infected after: {}", afterCountInfected)
    logger.info("Recovered: {}", afterCountRecovered)
    logger.info("Susceptible: {}", afterCountSusceptible)
  }

  private def getSusceptibleCount(context: Context) = {
    context.graphProvider.fetchCount("Person", "infectionState" equ Susceptible)
  }

  private def getInfectedCount(context: Context): Int = {
    context.graphProvider.fetchCount("Person", ("infectionState" equ Infected))
  }

  private def getRemovedCount(context: Context) = {
    context.graphProvider.fetchCount("Person", "infectionState" equ Removed)
  }
}

The output message on running the code is

_images/OutputFile_msg.png

Expanding the Network

Ealier we had one location class which was the House. In this section we increase the location classes to House, Office, and School. Every person has a unique house and either an office or a school and this categorized on the basis of age.

Implementing multiple houses, offices, and schools

As mention while creating the House.scala class, we mentioned that each of the locations will require a separate class. In addition to the new location classes, the person class needs to updated to establish the relationships.

This scala class defines the relationship betweeen the agent of type Person and Office. Again since there are numerous offices, the datatype required is Long.

package sir

import com.bharatsim.engine.models.Network

case class Office(id: Long) extends Network {
  addRelation[Person]("EMPLOYER_OF")

  override def getContactProbability(): Double = 1
}

The main file doesnt need major alterations, but the changes that have to be implemented are crucial conceptually and for the program to give the correct output. The majority of the changes are in two areas which are

  • Categorization of people: We have different locations in the network but only one type of Person. We need to make a distinction and categorize the individuals to send them to different locations. In this section, the categorization is done on the basis of age; any over the age of 18 works in an office and anyone under the age of 18 goes to a school. After creating these different people, we need to define the relationship between the people and their respective nodes. All these changes are made in the csvDataExtractor.

Note

The age of the citizens are provided in the input csv file.

  • createSchedules: Now that we have defined office-goers and school-goers, we need to decide their schedules and timings.

The csvDataExtractor function is the same and changes are made after the nodes (house, citizen) and relationship (house and person) is defined. Regardless of the age of the individual, they still have a house that they are associated to and therefore no changes are required when defining the aforementioned nodes and relationships. The next part is adding new nodes and relationships for individuals and their additional network and this is rather straightforward. An if condition is used to categorize on the basis of age and in the conditional block the relationships and nodes are added, similar to the house and citizen case.

if (age >= 18) {
  val office = Office(officeId)
  val worksAt = Relation[Person, Office](citizenId, "WORKS_AT", officeId)
  val employerOf = Relation[Office, Person](officeId, "EMPLOYER_OF", citizenId)

  graphData.addNode(officeId, office)
  graphData.addRelations(worksAt, employerOf)
} else {
  val school = School(schoolId)
  val studiesAt = Relation[Person, School](citizenId, "STUDIES_AT", schoolId)
  val studentOf = Relation[School, Person](schoolId, "STUDENT_OF", citizenId)

  graphData.addNode(schoolId, school)
  graphData.addRelations(studiesAt, studentOf)
}

After this distinction has been made, the changes in schedules have to be made. Employee and student schedule are just when they leave for their the house and when they return. First we need to define an hour to be myTick and there are 24 hours in myDay. Before create24HourSchedules can be made, myTick and myDay needs to be defined outside the main function.

private val myTick: ScheduleUnit = new ScheduleUnit(1)
private val myDay: ScheduleUnit = new ScheduleUnit(myTick * 24)

With these values defined, create24HourSchedules can be made. However, when there are more than one schedules running, there needs to be a priority list that needs to be made. In this case, Student and Employee schedules are independent of each other so a either schedules can be prioritized over the other. In later cases, quarantine will be introduced where individuals will stay at their house the whole time and this gets priority over office and school schedules.

private def create24HourSchedules()(implicit context: Context): Unit = {
  val employeeSchedule = (myDay, myTick)
    .add[House](0, 8)
    .add[Office](9, 17)
    .add[House](18,23)

  val studentSchedule = (myDay, myTick)
    .add[House](0, 8)
    .add[Office](9, 16)
    .add[House](17, 23)

  registerSchedules(
    (employeeSchedule, (agent: Agent, _: Context) => agent.asInstanceOf[Person].age >= 18, 1),
    (studentSchedule, (agent: Agent, _: Context) => agent.asInstanceOf[Person].age < 18, 2)
  )
}

Note

The timings of departure and return are to be made in the 24 hour format.

The whole main file code is

package sir

import java.util.Date
import com.bharatsim.engine.ContextBuilder._
import com.bharatsim.engine._
import com.bharatsim.engine.actions.StopSimulation
import com.bharatsim.engine.basicConversions.decoders.DefaultDecoders._
import com.bharatsim.engine.basicConversions.encoders.DefaultEncoders._
import com.bharatsim.engine.dsl.SyntaxHelpers._
import com.bharatsim.engine.execution.Simulation
import com.bharatsim.engine.graph.ingestion.{GraphData, Relation}
import com.bharatsim.engine.graph.patternMatcher.MatchCondition._
import com.bharatsim.engine.listeners.{CsvOutputGenerator, SimulationListenerRegistry}
import com.bharatsim.engine.models.Agent
import com.bharatsim.engine.utils.Probability.biasedCoinToss
import com.bharatsim.examples.epidemiology.sir.InfectionStatus._
import com.typesafe.scalalogging.LazyLogging

object Main extends LazyLogging {
  private val initialInfectedFraction = 0.01

  private val myTick: ScheduleUnit = new ScheduleUnit(1)
  private val myDay: ScheduleUnit = new ScheduleUnit(myTick * 24)

  def main(args: Array[String]): Unit = {

    var beforeCount = 0
    val simulation = Simulation()

    simulation.ingestData(implicit context => {
      ingestCSVData("citizen10k.csv", csvDataExtractor)
      logger.debug("Ingestion done")
    })

    simulation.defineSimulation(implicit context => {
      create24HourSchedules()

      registerAction(
        StopSimulation,
        (c: Context) => {
          getInfectedCount(c) == 0
        }
      )

      beforeCount = getInfectedCount(context)

      registerAgent[Person]

      val currentTime = new Date().getTime

      SimulationListenerRegistry.register(
        new CsvOutputGenerator("src/main" + currentTime + ".csv", new SIROutputSpec(context))
      )
    })

    simulation.onCompleteSimulation { implicit context =>
      printStats(beforeCount)
      teardown()
    }

    val startTime = System.currentTimeMillis()
    simulation.run()
    val endTime = System.currentTimeMillis()
    logger.info("Total time: {} s", (endTime - startTime) / 1000)
  }

  private def create24HourSchedules()(implicit context: Context): Unit = {
    val employeeSchedule = (myDay, myTick)
      .add[House](0, 8)
      .add[Office](9, 17)
      .add[House](18,23)

    val studentSchedule = (myDay, myTick)
      .add[House](0, 8)
      .add[Office](9, 16)
      .add[House](17, 23)

    registerSchedules(
      (employeeSchedule, (agent: Agent, _: Context) => agent.asInstanceOf[Person].age >= 18, 1),
      (studentSchedule, (agent: Agent, _: Context) => agent.asInstanceOf[Person].age < 18, 2)
    )
  }

  private def csvDataExtractor(map: Map[String, String])(implicit context: Context): GraphData = {

    val citizenId = map("Agent_ID").toLong
    val age = map("Age").toInt
    val initialInfectionState = if (biasedCoinToss(initialInfectedFraction)) "Infected" else "Susceptible"

    val homeId = map("HHID").toLong
    val schoolId = map("school_id").toLong
    val officeId = map("WorkPlaceID").toLong

    val citizen: Person = Person(
      citizenId,
      age,
      InfectionStatus.withName(initialInfectionState),
      0
    )

    val home = House(homeId)
    val staysAt = Relation[Person, House](citizenId, "STAYS_AT", homeId)
    val memberOf = Relation[House, Person](homeId, "HOUSES", citizenId)

    val graphData = GraphData()
    graphData.addNode(citizenId, citizen)
    graphData.addNode(homeId, home)
    graphData.addRelations(staysAt, memberOf)

    if (age >= 18) {
      val office = Office(officeId)
      val worksAt = Relation[Person, Office](citizenId, "WORKS_AT", officeId)
      val employerOf = Relation[Office, Person](officeId, "EMPLOYER_OF", citizenId)

      graphData.addNode(officeId, office)
      graphData.addRelations(worksAt, employerOf)
    } else {
      val school = School(schoolId)
      val studiesAt = Relation[Person, School](citizenId, "STUDIES_AT", schoolId)
      val studentOf = Relation[School, Person](schoolId, "STUDENT_OF", citizenId)

      graphData.addNode(schoolId, school)
      graphData.addRelations(studiesAt, studentOf)
    }

    graphData
  }

  private def printStats(beforeCount: Int)(implicit context: Context): Unit = {
    val afterCountSusceptible = getSusceptibleCount(context)
    val afterCountInfected = getInfectedCount(context)
    val afterCountRecovered = getRemovedCount(context)

    logger.info("Infected before: {}", beforeCount)
    logger.info("Infected after: {}", afterCountInfected)
    logger.info("Recovered: {}", afterCountRecovered)
    logger.info("Susceptible: {}", afterCountSusceptible)
  }

  private def getSusceptibleCount(context: Context) = {
    context.graphProvider.fetchCount("Person", "infectionState" equ Susceptible)
  }

  private def getInfectedCount(context: Context): Int = {
    context.graphProvider.fetchCount("Person", ("infectionState" equ Infected))
  }

  private def getRemovedCount(context: Context) = {
    context.graphProvider.fetchCount("Person", "infectionState" equ Removed)
  }
}

Introduction of Social Interventions

BharatSim allows the user the possibility to implement a variety of social interventions that modify disease spread at the level of individuals, like quarantines, lockdowns, or vaccination drives. In this section, we will describe a simple implementation of a quarantine.

Quarantine

Quarantine can be brought into effect by forcing a schedule onto the people where everyone stays at their respective house. In create24HourSchedules everyone can be made to stay at home from 0 to 23, and this can be given the number 1 priority. When brought into effect, the school and office schedules will be ignored and the quarantine schedules will be abided by.

private def create24HourSchedules()(implicit context: Context): Unit = {
  val employeeSchedule = (myDay, myTick)
    .add[House](0, 8)
    .add[Office](9, 17)
    .add[House](18,23)

  val studentSchedule = (myDay, myTick)
    .add[House](0, 8)
    .add[Office](9, 16)
    .add[House](17, 23)

  val quarantinedSchedule = (myDay, myTick)
    .add[House](0, 23)

  registerSchedules(
    (quarantinedSchedule, (agent: Agent, _: Context) => agent.asInstanceOf[Person].isInfected, 1),
    (employeeSchedule, (agent: Agent, _: Context) => agent.asInstanceOf[Person].age >= 18, 2),
    (studentSchedule, (agent: Agent, _: Context) => agent.asInstanceOf[Person].age < 18, 3)
  )
}

FSM in SIR

Futher information can be found Finite State Machine

In this case of Finite State Machine, the abstract machine will be the agents in the population. Each of these agents can be either Susceptible, Infected or Recovered, and there is a well defined procedure on moving from one state to another. When the FSM is implemented, the Transition condition will be dictated by the probability of contracting the infection or recovering from the infection. In more complicated systems, there can be 2 or more transitions are possible, for example from Exposed to Asymptomatic or Symptomatic states. These transitions will be dictated by there respective probabilities and again only one of these transitions can take place.

Earlier we had introduced Disease Dynamics in the form of behaviours, and these dictated whether an agent would be isInfected or isRecovered. As discussed earlier, in the FSM Transition will dictate whether an agent is in InfectedState or RecoveredState.

To introduce a Finite State Machine, we need to make the following changes:

  1. Define disease states, and define transitions between them

  2. Modify our agents to now be extensions of the StatefulAgent class, instead of the Agent class like before.

package sir

object Disease {
  final val beta: Double = 0.3
  final val lastDay: Int = 12
  final val lamda: Double = 0.14
  final val dt: Double = 0.5
}

Right click on the folder, and hover over the Refactor option, and then click on copy classes. Rename these sets of classes as FSMsir.

Error

If the name appears as FSMsir.sir then simply rename the file through refactor as FSMsir

The Person extends the Agent class but now that we are re-defining how a person is thought off, we need to extend a pre-defined class called StatefulAgent. There is no need to import another package, it was added in the code above. Create a new package in the current and name it DiseaseStates, and create case classes called SusceptibleState, InfectedState, and RecoveredState.

In the DiseaseStates classes, import the following packages,

import com.bharatsim.engine.Context
import com.bharatsim.engine.basicConversions.decoders.DefaultDecoders._
import com.bharatsim.engine.basicConversions.encoders.DefaultEncoders._
import com.bharatsim.engine.fsm.State
import com.bharatsim.engine.graph.patternMatcher.MatchCondition._
import com.bharatsim.engine.models.{Network, StatefulAgent}
import com.bharatsim.engine.utils.Probability.biasedCoinToss
import FSMsir.InfectionStatus._
import FSMsir.{Disease, Person}

For each of the classes also extend the State Class. What we aim to achieve in these classes is have a defination of what it means to be of that State and a Transition out of that State.

Note

It is important to note that we define the probability to leave that State and not enter the State. That will be defined in the previous State to the current.

By doing so, we can remove a major portion of the code written in the Person class, since that was the governing the disease dynamics. It is more convenient to start by defining the Transition. The syntax is as follows,

addTransition(
when = ,
to = context =>
)

addTransition requires two parameters, when to execute the Transition and where does the agent go. The former is Boolean while the latter is a State. To tackle the when parameter we can define a function called shouldBeInfected, which does the same thing as checkforInfection in the Person class. As to where the agent will go after the Transition, that is the InfectedState we have just written. The Transition will be the following,

addTransition(
when = shouldBeInfected,
to = context => InfectedState()
)

Now it comes to defining the shouldbeInfected function, and this can be done by updating the checkforInfection function, however I use a different approach. This has incorporated PerTickCache method to reduce computational time. I will briefly explain the advantages of this method of computation. More often that not, there are multiple agents present at one location at any given tick and the current simulation calculates quantities like infectedCount, infectedNeighbourCount for each and every of these agents. At every Tick, the system has become static and the information of the location does not change, and it is becomes tedious to calculate all these quantities over and over again. PerTickCache calculates the information about the location once, and stores the information. If another agent belongs to the locations whose information was previously computed, then the stored information is utilized and if there is no information present, then it calculates and stores it for any other agent who might be present here. After the Tick has been completed, then it deletes the information. If there are N locations, then there will be a maximum of N times these quantities will be calculated.

def shouldBeInfected(context: Context, agent: StatefulAgent): Boolean = {
  if (agent.activeState == SusceptibleState()) {
    val infectionRate = Disease.beta
    val dt = Disease.dt

    val schedule = context.fetchScheduleFor(agent).get

    val currentStep = context.getCurrentStep
    val placeType: String = schedule.getForStep(currentStep)

    val places = agent.getConnections(agent.getRelation(placeType).get).toList
    if (places.nonEmpty) {
      val place = places.head
      val decodedPlace = agent.asInstanceOf[Person].decodeNode(placeType, place)

      val infectedFraction = fetchInfectedFraction(decodedPlace, placeType, context)
      return biasedCoinToss(infectionRate * infectedFraction * dt)
    }
  }
  false
}

This function is every similar to checkforInfection except for the conversion from Agent to StatefulAgent. Here the infectedFraction is not calculated, instead a value from another function is obtained. This is where PerTickCache is implemented.

private def fetchInfectedFraction(decodedPlace: Network, place: String, context: Context): Double = {
  val cache = context.perTickCache

  val tuple = (place, decodedPlace.internalId)
  cache.getOrUpdate(tuple, () => fetchFromStore(decodedPlace)).asInstanceOf[Double]
}

The above is a hashmap, which requires a tuple as a key which is unique for every location. The key which is a tuple that stores the place and the internalId of the place. This is fed in getOrUpdate, which looks into the stored memory to see if any information about the place can be found. If there exist some prior information, then it gets the information. If there is no prior information, then it calculates the values and updates it so the next time it will not have to calculate. The symbol () means that there is no information is present, and the computer is asked to use the function fetchFromStore to find the infected number. This is the same code as the one in Person class.

private def fetchFromStore(decodedPlace: Network): Double = {
  val infectedPattern =
    ("infectionState" equ Infected)
  val total = decodedPlace.getConnectionCount(decodedPlace.getRelation[Person]().get)

  total
}

These are the things that need to added to SusceptibleState class. From the agent Transitions to InfectedState. Again it is easier to add the Transition first.

addTransition(
  when = checkForRecovered,
  to = context => RecoveredState()
)

The function checkForRecovered is just a biasedCoinToss with the appropriate probabilities.

def checkForRecovered(context: Context, agent: StatefulAgent): Boolean = {
  return biasedCoinToss(Disease.lamda * Disease.dt)
}

This is all for InfectedState. Nothing needs to be added for RecoveredState since they cant participate in the dynamics or Transition out of the State. However, if we were to model a system with waning immunity to a disease - for example, an “SIRS” model where recovered individuals transition back to the Susceptible state – we will need to include the dynamics of this in the RecoveredState.

Miscellaneous

Assembling an Executable .jar file

A jar file is an executable java file that can be used to run your code. Instead of importing and building the whole project, it is more efficient to convert the source code into a jar and import it as a library. This way the environment will be clutter free and all the classes and functions can be easily imported.

In order to assemble the code into an executable jar that can be used as a library, open the project in IntelliJ and go to the sbt shell, and type assembly.

_images/jar-doc-1.png

Alternatively, you could also navigate to the root folder of your project in a standard UNIX terminal and type sbt assembly. Once this command has successfully run, a new jar file will appear in the target/scala-2.##/ within the root folder of your project.

Importing a .jar file as a Library

As mentioned in the above section, importing the whole source code as a jar file is highly efficient and recommended, especially for a novice coder. In the root folder of the project, create a new folder called lib unless one already exists. Copy the jar file and paste it into the lib folder.

  1. Go to file option in IntelliJ and click on the option “Project structure”.

  2. Select “modules” on the left side panel, and select the “dependencies” option. A list of files should appear.

  3. There should be a “+” option and select the Library followed by Java.

  4. Select the required jar file in the lib folder. Click on Apply and then ok button to import the jar file.

Error

If the package or a class from the library is not available, then make sure the jar file has been imported successfully. Once common source of error is not clicking on the Apply button before selecting ok as mentioned in Step 4.

Importing packages and classes can be done as before using the same syntax.

Using args in main method

While using the run command on the sbt shell, one can pass in some string arguments. These arguements can be called by the main function.

For example, one might want to change the name of the output file everytime they run the code. Instead of changing the output name manually each time, one can write a code like the one given below.

def main(args: Array[String]): Unit = {
  outputName = args(0)
  }

One could use the outputName and use it in the name of the csv file, for eg.

SimulationListenerRegistry.register(
    new CsvOutputGenerator(outputName+".csv", new SEIROutputSpec(context))
  )

To implement this, one must go to the sbt shell and type run "my_sir_model". If this were to be used on the code-block above, the output csv file will be named as my_sir_model.csv. However, if one runs the file without specifiying the arguement, it will show an error:

_images/jar-doc-5.png

One can also run a Main file by creating a .jar file, as described above and then running java -jar file.jar [ arguments ]

Saving location-level information from the simulation

Tip

Before reading this section, it’s recommended that you read the basics in the Outputs section.

In this section, we discuss how to create the following output:

  • We require the function to accept a workplace type as a string (i.e. “Home”, “Office”, “School” in our model), and create a CSV file with the output.

  • Each file should have the Location type (“Home”, etc.) the location ID, and the number of people per location broken up by age-group (<18, 18-45, 45-60, 60+)

Let’s call our user-defined function myCsvOutputSpec for now. First, we’ll create a scala class for it which is an extention of the CsvSpecs trait:

class myCsvOutputSpec(placeType: String, context: Context) extends CSVSpecs {}

Next, we ovverride the getHeaders function with the appropriate list of headers:

override def getHeaders: List[String] = List("PlaceType", "LocationID", "N <18", "N 18-45", "N 45-60", "N >60")

Before overriding getRows, let’s write down two functions which we’ll be needing inside of it. The first is the decodeNode method, which converts a GraphNode to a Node

def decodeNode(classType: String, node: GraphNode): Node = {
  classType match {
    case "House" => node.as[House]
    case "Office" => node.as[Office]
    case "School" => node.as[School]
  }
}

The next method is called getId, and it retrieves the location ID of a GraphNode. In our program, the House, Office and School classes all have an attribute called id, so this function is designed to return that attribute.

def getId(classType: String, node: GraphNode) : Long = {
  classType match {
    case "House" => node.as[House].id
    case "Office" => node.as[Office].id
    case "School" => node.as[School].id
  }
}

Caution

Looking at this function, you may think it’s unnecessary: it looks almost identical to decodeNode! Why not just use decodeNode(classType, node).id? In that case, however, note that decodeNode returns a Node, which does not have an id attribute.

By playing around with the function, you may find out that the GraphNode attribute does have an id: so why not just write the function to return node.id? The GraphNode.id attribute is a completely different number from the location ID, which is used to identify the node on the graph. As such, while the code will compile and run, the output under LocationID will have different results from what you’d expect.

Now, we can start to write down our getRows method. We want to be able to initialize a large list, every component of which is a list containing a row of the CSV file. While it sounds tempting to first initialize an empty list, and add lists to it one at a time, that is not possible in scala. This is because the List datatype is immutable - although you can define a list just fine, it cannot be changed after. We can get around this by using the ListBuffer datatype, which has a lot of useful methods.

override def getRows(): List[List[Any]] = {

  val rows = ListBuffer.empty[List[String]]

}

Next we get all the nodes of the correct placeType (which, remember, was a string that the function accepts as an argument)

val locations = context.graphProvider.fetchNodes(placeType)

Iterating over each location, which we call oneLocation:

locations.foreach(oneLocation => {})

We generate a decodedLoc and locId using our decodeNode and getId functions respectively

val decodedLoc = decodeNode(placeType, oneLocation)
val locId = getId(placeType, oneLocation).toString

Note

We convert locId to a string, as it’s what we need to fill out as the second element of the row.

We then calculate the number of people in each age group who are associated with the location: This is done with getConnectionCount, where we feed in the relation between the location and the person, and then the age-requirement. We then convert the numbers to strings.

val N_0_18 = decodedLoc.getConnectionCount(decodedLoc.getRelation[Person]().get,
  "age" lt 18).toString
val N_18_45 = decodedLoc.getConnectionCount(decodedLoc.getRelation[Person]().get,
  ("age" gte 18) and ("age" lt 45)).toString
val N_45_60 = decodedLoc.getConnectionCount(decodedLoc.getRelation[Person]().get,
  ("age" gte 45) and ("age" lt 60)).toString
val N_60_100 = decodedLoc.getConnectionCount(decodedLoc.getRelation[Person]().get,
  "age" gte 60).toString

Now, we add this row to rows, the ListBuffer object

rows.addOne(List(placeType, locId, N_0_18, N_18_45, N_45_60, N_60_100))

Finally, outside of the iterator, we convert the ListBuffer to a List and return it

rows.toList

Putting it all together, the class is

class myCsvOutputSpec(placeType: String, context: Context) extends CSVSpecs {

    override def getHeaders: List[String] = List("PlaceType", "LocationID", "N_<18", "N_18-45", "N_45-60", "N_>60")


    override def getRows(): List[List[Any]] = {

        val rows = ListBuffer.empty[List[String]]

        val locations = context.graphProvider.fetchNodes(placeType)

        locations.foreach(oneLocation => {
        val decodedLoc = decodeNode(placeType, oneLocation)
        val locId = getId(placeType, oneLocation).toString
        val N_0_18 = decodedLoc.getConnectionCount(decodedLoc.getRelation[Person]().get, "age" lt 18).toString
        val N_18_45 = decodedLoc.getConnectionCount(decodedLoc.getRelation[Person]().get, ("age" gte 18) and ("age" lt 45)).toString
        val N_45_60 = decodedLoc.getConnectionCount(decodedLoc.getRelation[Person]().get, ("age" gte 45) and ("age" lt 60)).toString
        val N_60_100 = decodedLoc.getConnectionCount(decodedLoc.getRelation[Person]().get, "age" gte 60).toString

        rows.addOne(List(placeType, locId, N_0_18, N_18_45, N_45_60, N_60_100))
        })
        rows.toList
    }

    def decodeNode(classType: String, node: GraphNode): Node = {
        classType match {
        case "House" => node.as[House]
        case "Office" => node.as[Office]
        case "School" => node.as[School]
        }
    }

    def getId(classType: String, node: GraphNode) : Long = {
        classType match {
        case "House" => node.as[House].id
        case "Office" => node.as[Office].id
        case "School" => node.as[School].id
        }
    }

}

Tip

If you want to use this code snippet, be sure to import the following

import com.bharatsim.engine.Context
import com.bharatsim.engine.basicConversions.decoders.DefaultDecoders._
import com.bharatsim.engine.basicConversions.encoders.DefaultEncoders._
import com.bharatsim.engine.graph.GraphNode
import com.bharatsim.engine.graph.patternMatcher.MatchCondition._
import com.bharatsim.engine.listeners.CSVSpecs
import com.bharatsim.engine.models.Node
import scala.collection.mutable.ListBuffer

As we only need to call this function once after data ingestion, we add the following inside simulation.defineSimulation:

var outputGenerator = new CsvOutputGenerator("output.csv", new myCsvOutputSpec("House", context))
outputGenerator.onSimulationStart(context)
outputGenerator.onStepStart(context)
outputGenerator.onSimulationEnd(context)

The output should be of the form

output.csv

PlaceType

LocationID

N <18

N 18-45

N 45-60

N >60

House

1956

0

0

3

3

House

1762

0

3

1

2

House

680

1

1

0

1

House

1584

1

1

0

0

House

1730

0

4

0

0

House

2096

1

2

0

0

House

1903

0

2

0

1

House

1414

0

0

1

0

House

1087

1

0

2

2

House

652

1

0

0

2

Other examples

Caution

Before having a look at these examples, it’s highly recommended that you go through the Writing your First Program section!

Houses on a lattice

The model discussed earlier has no sense of geography: each house and office is practically isolated from all the others, as you can only infect people in your location. A simple way to try to introduce this into the model is by using houses on a grid, with the following characteristics:

  • Houses are arranged periodically in rows and columns

  • Each house has 4 nearest neighbours - namely the houses north, south, east, and west of them

  • Assume that people in a house can only infect their nearest neighbours

  • The inhabitants of a house are less likely to infect their neighbours than they are to their housemates

The house arrangement

The lattice of houses - each square represents a house. Here, houses 2,6,8 and 12 are neighbours of 7, and could potentially infect or get infected by it.

To simplify what happens at the boundaries, we implement periodic boundary conditions. This means the left edge and the right edge are identified, and likewise for the top and bottom edges of the lattice. Apart from being easy to implement, a system with periodic bondary conditions can be thought of as one which feels the influence of it’s environment, rather than being isolated.

Depiction of periodic boundary conditions

With periodic boundary conditions, house 24 is now to the left of 20, and 0 is below 20.

We use a 50 × 50 lattice, which totals to 2500 houses. We can also add an intervention to the model, to make it more realistic. We divide the lattice into larger ‘blocks’ of 100 houses each, so that a 5 × 5 arrangement of 25 blocks gives us our original lattice.

Note

The implementation below uses a block attribute in the house class. The blocks are numbered in the same way the houses are.

We can keep track of the number of infected in each block, and block it off if the infected count rises too high.

Depiction of a block being blockaded

If we split a 4 × 4 lattice into 4 blocks, we can see the boundary cordoned off if block 1 has a large number of infected.

We shall cordon off a block in the following way:

  • If infectedCountToTriggerBlockade is above a threshold (we use 30), then the block needs to be closed off

  • The block is closed off for blockadeDuration from the instant where the threshold is breached, regardless of the current levels of infected

  • After the block comes out of the blockade, it cannot get locked down again for blockadeCooldown days, regardless of the current levels of infected

Rather than having the initial infected people spread around, we’ll set them to be clustered in an area to better see the spread. We’ll set the four corners as the infected area: despite looking like 4 separate regions they’re actually all connected due to the periodic boundary conditions.

The model is built over an SIR model, with most of the changes being made to three classes: Main, SusceptibleState and GISOutputSpec.

There are several changes made to this class, the most important one being the addition of the intervention.

We first initialize isBlockadedList, a boolean array of length 25. If isBlockadedList[i] == True, that means that block number i is currently blockaded.

We use an IntervalBasedIntervention to set up our required intervention, and check for whether a block needs to be cordoned off every tick with whenActiveActionFunc.

Note

We start the intervention at tick 1, and end it at tick 5000. The simulation is assumed to have ended by then, as after numerous trials it never seemed to reach 5000 ticks. However, this is not good practice: we ideally should force the simulation to end when the intervention does, by adding the required condition in StopSimulation.

When ingesting the input data, we define a NEIGHBOURS relation between two neighbouring houses. However, we cannot register a relation between two nodes if only one of them has been registered. Thus, for every row in the input csv file, we register not just the house mentioned on the row but also all of it’s neighbours.

Hint

If you refer back to the final note in the section on Reading inputs from a synthetic population, you’ll see that registering an identical node multiple times does not lead to the node being duplicated in the graph.

What does what?

Here’s a breakdown of the primary methods in the Main class, and what they do.

  • blockadeBlock: Defines and registers the intervention described above

  • create12HourSchedules: Defines and registers the agent registerSchedules

  • csvDataExtractor: Creates the graph using the CSV input file

  • getLeftNeighbour, getRightNeighbour, getUpNeighbour, getDownNeighbour: Take a houses’ ID and returns the ID of the appropriate neighbour, after taking periodic boundary conditions into account

  • getHouseBlock: Returns the block number a given house is part of

The source code for these classes can be found below:

package com.bharatsim.examples.epidemiology.latticeHousesModel

import com.bharatsim.engine.ContextBuilder._
import com.bharatsim.engine._
import com.bharatsim.engine.actions.StopSimulation
import com.bharatsim.engine.basicConversions.decoders.DefaultDecoders._
import com.bharatsim.engine.basicConversions.encoders.DefaultEncoders._
import com.bharatsim.engine.dsl.SyntaxHelpers._
import com.bharatsim.engine.execution.Simulation
import com.bharatsim.engine.graph.ingestion.{GraphData, Relation}
import com.bharatsim.engine.graph.patternMatcher.MatchCondition._
import com.bharatsim.engine.intervention.{IntervalBasedIntervention, SingleInvocationIntervention}
import com.bharatsim.engine.listeners.{CsvOutputGenerator, SimulationListenerRegistry}
import com.bharatsim.engine.models.{Agent, Node}
import com.bharatsim.examples.epidemiology.latticeHousesModel.DiseaseStates.{InfectedState, SusceptibleState}
import com.bharatsim.examples.epidemiology.latticeHousesModel.InfectionStatus._
import com.typesafe.scalalogging.LazyLogging

import java.util.Date

object Main extends LazyLogging {

  final val numberOfTicksInADay: Int = 2
  final val dt: Double = 1/numberOfTicksInADay.toFloat

  private val myTick: ScheduleUnit = new ScheduleUnit(1)
  private val myDay: ScheduleUnit = new ScheduleUnit(myTick * numberOfTicksInADay)

  var isBlockadedList = new Array[Boolean](25)

  def main(args: Array[String]): Unit = {
    var beforeCount = 0
    val simulation = Simulation()

    simulation.ingestData(implicit context => {
      ingestCSVData("citizen10kLattice.csv", csvDataExtractor)
      logger.debug("Ingestion done")
    })

    simulation.defineSimulation(implicit context => {
      create12HourSchedules()

      blockadeBlock

      registerAction(
        StopSimulation,
        (c: Context) => {
          getInfectedCount(c) == 0
        }
      )

      beforeCount = getInfectedCount(context)

      registerAgent[Person]

      val currentTime = new Date().getTime

      SimulationListenerRegistry.register(
        new CsvOutputGenerator("src/main/resources/GISInfectedoutput_"+currentTime+".csv", new GISOutputSpec(context))
      )
    })

    simulation.onCompleteSimulation { implicit context =>
      printStats(beforeCount)
      teardown()
    }

    val startTime = System.currentTimeMillis()
    simulation.run()
    val endTime = System.currentTimeMillis()
    logger.info("Total time: {} s", (endTime - startTime) / 1000)
  }

  private def blockadeBlock(implicit context: Context): Unit = {

    val interventionName = "blockade"
    val infectedCountToTriggerBlockade = 30
    val blockadeDuration = 7 * numberOfTicksInADay
    val blockadeCooldown = 7 * numberOfTicksInADay
    var ticksSinceBlockade = Array.fill(25){0}

    def perTickAction(context: Context): Unit = {
      for (i <- 0 to 24)  {

        if (ticksSinceBlockade(i) == blockadeDuration) {
          isBlockadedList(i) = false
        }

        if (ticksSinceBlockade(i) >= blockadeDuration + blockadeCooldown) {
          var infectedCountPerBlock: Long = 0
          var nodesInBlock = context.graphProvider.fetchNodes("House", "block" equ i)
          nodesInBlock.foreach(blockNode => {
            var tempvariable = fetchInfectedAndTotalPerLocation(blockNode.as[House], "House", context)
            infectedCountPerBlock += tempvariable._1.toLong
          }
          )

          if (infectedCountPerBlock >= infectedCountToTriggerBlockade) {
            isBlockadedList(i) = true
            ticksSinceBlockade(i) = 0
          }

        }
        else ticksSinceBlockade(i) += 1
      }
    }

    def fetchInfectedAndTotalPerLocation(node: Node, placeType: String, context: Context): (Double, Double) = {
      val cache = context.perTickCache
      val uniquekey = (placeType, node.internalId)
      cache.getOrUpdate(uniquekey, () => computeInfectedAndTotalPerLocation(node)).asInstanceOf[(Double, Double)]
    }

    def computeInfectedAndTotalPerLocation(node: Node): (Double, Double) = {
      val totalNeighbourCount = node.getConnectionCount(node.getRelation[Person]().get)
      if (totalNeighbourCount == 0)
        return (0d, 1)  // toDo change to (0,0), add check for dividing by 0
      val infectedNeighbourCount = node.getConnectionCount(node.getRelation[Person]().get,
        "infectionState" equ Infected)
      return (infectedNeighbourCount.toDouble, totalNeighbourCount.toDouble)
    }

    val intervention =
      IntervalBasedIntervention(interventionName, 1, 5000, whenActiveActionFunc = perTickAction)

    registerIntervention(intervention)
  }

  private def create12HourSchedules()(implicit context: Context): Unit = {

    val stayHomeSchedule = (myDay, myTick)
      .add[House](0, 1)

    registerSchedules(
      (stayHomeSchedule, (agent: Agent, _:Context) => agent.asInstanceOf[Person].age > 0, 1)
    )
  }

  private def csvDataExtractor(map: Map[String, String])(implicit context: Context): GraphData = {

    val citizenId = map("Agent_ID").toLong
    val age = map("Age").toInt

    val homeId = map("HHID").toLong
    val schoolId = map("school_id").toLong
    val officeId = map("WorkPlaceID").toLong
    val houseLatitude = map("H_Lat").toString
    val houseLongitude = map("H_Lon").toString

    val initialInfectionState = if ((houseLatitude=="0" || houseLatitude=="1" || houseLatitude=="49" ||
      houseLatitude=="2" || houseLatitude=="48") && (houseLongitude=="0" || houseLongitude=="1" || houseLongitude=="49" ||
      houseLongitude=="2" || houseLongitude=="48")) "Infected" else "Susceptible"

    val citizen: Person = Person(
      citizenId,
      age,
      houseLatitude,
      houseLongitude,
      InfectionStatus.withName(initialInfectionState),
      0,
      getInitialRecoveryTick(initialInfectionState)
    )

    if (initialInfectionState == "Susceptible") {
      citizen.setInitialState(SusceptibleState())
    }
    else
      citizen.setInitialState(InfectedState())

    val home = House(homeId, getHouseBlock(homeId))
    val staysAt = Relation[Person, House](citizenId, "STAYS_AT", homeId)
    val memberOf = Relation[House, Person](homeId, "HOUSES", citizenId)

    val neighboursLeft = Relation[House, House](homeId, "NEIGHBOURS", getLeftNeighbour(homeId))
    val neighboursRight = Relation[House, House](homeId, "NEIGHBOURS", getRightNeighbour(homeId))
    val neighboursUp = Relation[House, House](homeId, "NEIGHBOURS", getUpNeighbour(homeId))
    val neighboursDown = Relation[House, House](homeId, "NEIGHBOURS", getDownNeighbour(homeId))

    val graphData = GraphData()
    graphData.addNode(citizenId, citizen)
    graphData.addNode(homeId, home)
    graphData.addRelations(staysAt, memberOf)

    var lHomeId = getLeftNeighbour(homeId)
    var rHomeId = getRightNeighbour(homeId)
    var uHomeId = getUpNeighbour(homeId)
    var dHomeId = getDownNeighbour(homeId)

    graphData.addNode(lHomeId, House(lHomeId, getHouseBlock(lHomeId)))
    graphData.addNode(rHomeId, House(rHomeId, getHouseBlock(rHomeId)))
    graphData.addNode(uHomeId, House(uHomeId, getHouseBlock(uHomeId)))
    graphData.addNode(dHomeId, House(dHomeId, getHouseBlock(dHomeId)))

    graphData.addRelations(staysAt, memberOf)
    graphData.addRelations(neighboursLeft, neighboursRight, neighboursUp, neighboursDown)

    if (age >= 25) {
      val office = Office(officeId)
      val worksAt = Relation[Person, Office](citizenId, "WORKS_AT", officeId)
      val employerOf = Relation[Office, Person](officeId, "EMPLOYER_OF", citizenId)

      graphData.addNode(officeId, office)
      graphData.addRelations(worksAt, employerOf)
    } else {
      val school = School(schoolId)
      val studiesAt = Relation[Person, School](citizenId, "STUDIES_AT", schoolId)
      val studentOf = Relation[School, Person](schoolId, "STUDENT_OF", citizenId)

      graphData.addNode(schoolId, school)
      graphData.addRelations(studiesAt, studentOf)
    }

    graphData
  }

  private def getLeftNeighbour(houseID: Long) : Long = {
    if ((houseID + 1) % 50 == 0 ) {
      houseID + 1 - 50
    }
    else houseID + 1
  }

  private def getRightNeighbour(houseID: Long) : Long = {
    if (houseID % 50 == 0 ) {
      houseID - 1 + 50
    }
    else houseID - 1
  }

  private def getUpNeighbour(houseID: Long) : Long = {
    (houseID + 50) % 2500
  }

  private def getDownNeighbour(houseID: Long) : Long = {
    (houseID - 50 + 2500) % 2500
  }

  def getHouseBlock(houseID: Long) : Int = {
    val block_Lat = (houseID % 50) / 10
    val block_Lon = (houseID / 50) / 10
    5*block_Lat.toInt + block_Lon.toInt
  }


  private def printStats(beforeCount: Int)(implicit context: Context): Unit = {
    val afterCountSusceptible = getSusceptibleCount(context)
    val afterCountInfected = getInfectedCount(context)
    val afterCountRecovered = getRemovedCount(context)

    logger.info("Infected before: {}", beforeCount)
    logger.info("Infected after: {}", afterCountInfected)
    logger.info("Susceptible: {}", afterCountSusceptible)
    logger.info("Recovered: {}", afterCountRecovered)
  }

  private def getInitialRecoveryTick(state: String): Double = {
    if (state == "Susceptible") {
      0
    }
    else {
      numberOfTicksInADay*Disease.infectionDurationPDF.sample()
    }
  }

  private def getSusceptibleCount(context: Context) = {
    context.graphProvider.fetchCount("Person", "infectionState" equ Susceptible)
  }

  private def getInfectedCount(context: Context) = {
    context.graphProvider.fetchCount("Person", "infectionState" equ Infected)
  }

  private def getRemovedCount(context: Context) = {
    context.graphProvider.fetchCount("Person", "infectionState" equ Removed)
  }
}

Known Issues

java.util.NoSuchElementException on fetchActiveState

This issue will lead to the following output:

[ticks-loop-akka.actor.default-dispatcher-46] ERROR akka.actor.LocalActorRefProvider - guardian failed, shutting down system

Tip

Note that the 46 in the above error could be any number, as this error arises from parallelism

Below this, you will find one of the following errors:

java.util.NoSuchElementException: null
      at scala.collection.concurrent.TrieMap.apply(TrieMap.scala:878)
java.lang.Exception: Something went wrong, each StatefulAgent must have one active state all the times
      at com.bharatsim.engine.models.StatefulAgent.fetchActiveState(StatefulAgent.scala:69)

The error is not easily reproducible, and appears to occur randomly.

A preliminary investigation suggests a problem with the underlying TrieMap data structure. Edges on the map do not appear to be created, and cannot be referenced when called. Note that this behaviour only occurs if bharatsim is parallelized.

A potential workaround is to run bharatsim without parallelism. To do so, you need to edit the application.conf file, which can be found in /src/main/resources/. The first couple of lines should be

bharatsim {
  engine {
      execution {
          mode = "actor-based"
          mode = ${?EXECUTION_MODE}
          simulation-steps = 4000
          simulation-steps = ${?SIMULATION_STEPS}
          actor-based {
              num-processing-actors = 100
          }
      }

Line 4 needs to be changed to

mode = "no-parallelism"

This will fix the issue, but be aware that removing parallelism could lead to your code taking longer to run.

Note

For the latest on the error, please check the github issue