A simple virus contagion simulator

In March 2020, the Washington Post published an online paper that showed, through different simulations, how the spread of a virus evolves under different circumstances --such as the now well-known social distancing. Each of the simulations was made of a rectangular area, in the user's browser, where small circles could move around, bump into (and infect) each other, and eventually recover.

Depending on how many of them could move, the total and peak number of infected circles would follow different curves, such as the ones shown below.

This project is an attempt at reproducing this simulator, by using a variety of software tools developed at Laboratoire d'informatique formelle, a research lab in Computer Science from Université du Québec à Chicoutimi in Québec, Canada. Most notably, the project showcases two important libraries developed at LIF:

• The BeepBeep event stream processing engine
• The Synthia data structure generator

The rest of this Readme explains how the simulator has been built.

See a video of the simulator running.

Premise

The Post's simulator works as follows:

1. Circles are randomly placed in a rectangular "arena" and are initially moving in randomly selected directions. A varying number of circles can be made fixed: they don't move for the entire duration of the simulation.
2. A single circle is initially marked as "infected" (brown); all the other circles are "healthy" (light blue).
3. When two circles collide, they bounce off each other; moreover, if one is infected and the other is healthy, the healthy one becomes infected.
4. After a fixed amount of time, an infected circle becomes "recovered" (purple); it can neither become infected again, nor make other circles infected.

Simulating the "arena"

In our project, the management of the physics of colliding circles is done by the package virussim.physics. Inside this package, the Ball class implements a simple two-dimensional "ball", with a position and a speed vector, that can compute elastic collisions with other balls. The Arena is just a set of balls; its method update acts like a "tick" that moves the simulation one step forward: it goes through all Ball objects, updates their state, and determines if any collision occurs.

We shall not elaborate much on this part of the program, as most of the code in this package is borrowed and adapted from an answer on StackOverflow. One simple adaptation is that the arena can be asked to output its current state, in the form of a Map linking the unique ID of each ball with the corresponding ball instance.

The virussim package defines the Patient class, which is a descendent of Ball that can catch a virus when in contact with another infected ball. To this end, Patient has a member field keeping its Health state. More on that later.

Generating the initial state

The initial state of the arena is created by making multiple random "choices" for each patient. In order to generate this state, we make extensive use of the Picker interface provided by the Synthia library. A "picker" is any object that implements a method called pick(), which, when called, returns an object of a certain type. For example, RandomInteger is a picker that returns a randomly selected integer number on every call to its pick() method; the same for RandomFloat, RandomBoolean, and so on.

For example, to generate integers in the interval [5,15], a picker can be created in the following way:

f = new RandomInteger(5, 15);
System.out.println(f.pick());
System.out.println(f.pick());
...

The repeated calls to f.pick() print the sequence of integers generated by this picker, e.g.: 8, 6, 9, ...

Initial position

However, Synthia can generate objects that go beyond random scalars. To generate each patient's initial position, we use the PrismPicker: given two pickers p1 and p2, this picker generates a 2D vector where the first coordinate is given by asking p1, and the second coordinate is given by asking p2.

Depending on how these two values are picked, this can correspond to a different placement of the circles in the arena. Our program provides two options:

• In the first case, p1 and p2 perform an affine transformation of a uniformly distributed float in the interval [0,1]; in such a case, the circles are scattered uniformaly over the arena.
• In the second case, p1 and p2 perform a different affine transformation over a float that follows a Gaussian distribution. This causes the circles to be clustered towards the center of the arena.

For example, the following piece of code creates a picker for floats following a N(0;1) Gaussian distribution, and creates another picker that turns it into a N(H/2;H/6) distribution (for some value H):

gf = new GaussianFloat();
at = new AffineTransform.AffineTransformFloat(gh, H/6, H/2);

Initial velocity

Similarly, the initial direction of each patient is given by another 2D vector specifying its velocity. However, for the sake of elegance, it would be desirable that the speed of each circle be the same --that is, their velocity vector can have an arbitrary orientation, but must have the same modulus. This can be done using Synthia's HyperspherePicker, which does exactly that: it takes as input a vector length, and a picker that provides an arbitrary float value. This value is interpreted as an angle (in radians) in a polar coordinate system, thus generating vectors of a fixed modulus but with a potentially varying angle.

float modulus = 5f;
Picker<Float> pf = ...
HypershperePicker hp = new HyperspherePicker(modulus, pf);

Fixed or moving?

Another part of the patient's state that must be specified is whether its associated circle is fixed or moving. The Post's simulations test various proportions of moving circles. In our setup, the decision whether a patient is moving or not is provided by a straightforward RandomBoolean picker, which works like a biased coin toss; indeed, we can specify the probability that it returns true (meaning the patient is fixed) to whatever fraction we wish.

float probability = 0.3f;
RandomBoolean rb = new RandomBoolean(probability);

In this example, each call to rb.pick() will have a probability of 0.3 of returning true.

Updating health status

In the Post simulation, each sick patient transitions to the "recovered" state after a fixed amount of time (i.e. a predefined number of simulation steps). This could be handled with a simple counter variable inside the Patient class; but we chose to implement it using Synthia's picker interface to give us the possibility to explore other rules for recovery. Therefore, the Patient class has a member field m_healthPicker that contains an object implementing Picker<Health> --that is, a picker object that produces values of type Health. Each time this picker is queried, it is intended to return the patient's health state for the next simulation step.

First strategy: fixed recovery

The original recovery rule (recover after N steps) can be implemented using a special type of Synthia picker called Playback. This picker is instantiated with a list of values; it iterates through these values on each call to pick(), until it reaches the last one, which it repeatedly outputs from this point on.

In our case, we need to create a Playback picker that returns INFECTED N times, followed by RECOVERED indefinitely. It suffices to instantiate a Playback with the approrpiate values in its list:

Health[] h = new Health[N+1];
for (int i = 0; i < N; i++)
  h[i] = Health.INFECTED;
h[N] = Health.RECOVERED;
Playback pb = new Playback(h);

The Patient is programmed in such a way that, as long as its current state is INFECTED, its next state in the simulation is given by asking the picker. This indeed produces the desired behavior.

Second strategy: Markov chain

However, the fact that a patient uses a picker to update its health state means that this picker can be replaced by another one if we wish, without having to change anything in the Patient class. As a second strategy, we use another picker object provided by Synthia which is called the MarkovChain.

For the purpose of this project, a Markov chain can be seen as a finite-state machine where each transition from one state to the next is associated to a probability. In Synthia, the MarkovChain picker associates values to each state; repeated calls to its pick() method generate one possible run in the chain, where in each state, the next state is selected with a probability defined by the corresponding transition.

The Markov chain is a handy way of specifying the possible "behaviors" that something can have. We shall use it to define how a patient can transition from the INFECTED to the RECOVERED state. We will also take advantage of the fact that there can be multiple next states and add a second outcome for a sick patient, namely that it can transition to an additional health state we call DEAD (hopefully with a very small probability!). Graphically, this is illustrated by the following state machine:

In the infected (I) state, a patient has a probability pD of moving to the dead (D) state in the next simulation step, a probability pR of moving to the recovered state (R); the rest of the time, it remains in the infected state. We can also see that once a patient reaches states R or D, it then stays there forever. One can create this Markov chain in Synthia as follows:

Picker<Float> pf = ...
MarkovChain m = new MarkovChain(pf);
m.add(0, INFECTED);
m.add(1, RECOVERED);
m.add(2, DEAD);
m.add(0, 1, pR);
m.add(0, 2, pD);
m.add(0, 0, 1 - (pR + pD));
m.add(1, 1, 1);
m.add(2, 2, 1);

In the project, the HealthMarkovChain class simply creates such a Markov chain for given values of pR and pD. Note that the object also needs a Picker<Float>, which it uses to select the next transition to take on each call to pick().

Creating patients

Using such a model makes it possible to try the simulation with different recovery strategies, and different parameters for each strategy. It suffices to give a different Picker<Health> to the patient and restart the program. Each patient can even be given a different (randomly selected!) picker.

Equipped with these various pickers, it is now possible to generate patients. We do so by creating a PatientPicker, which is an object implementing Synthia's Picker<Patient> interface. In other words, it is a class that produces a new instance of Patient every time its pick() method is called. In order to do so, the PatientPicker requires:

• a Picker<float[]> to generate the position vector (e.g. the PrismPicker)
• a Picker<float[]> to generate the velocity vector (e.g. the HyperspherePicker)
• a Picker<Boolean> to determine if it is moving (e.g. the RandomBoolean)
• a Picker<Health> to manage its health status (e.g. the Playback or MarkovChain objects we just described)

Once instantiated in such a way, the PatientPicker is repeatedly called, and the patients it produces are added to the arena. The simulation is then ready to start.

Driving the simulation

So far, our simulation makes heavy use of objects provided by Synthia; but what about BeepBeep?

This event stream library will be used to drive the simulation and process its output in various ways. As you may know, BeepBeep is based on the concept of processor: a processor is a computing unit that receives data elements called events as its input, and produces other events as its output. In the present case, the initial source of events will be the arena --or rather the set of patients it contains. An event will be the Map object that associates each patient ID with the corresponding patient instance. A stream of events will be formed by the succession of such maps at each simulation step.

To this end, the Arena object must be turned into a special BeepBeep processor called a Source. This is the purpose of the ArenaSource class, which is a simple wrapper around an arena, so that its current state can be used as a source of BeepBeep events. Since ArenaSource is a BeepBeep processor, it can be connected to any other processor provided by the library.

Drawing the simulation

A first essential task of the simulator is to display the current position and state of each patient, as in the Post's "circle" simulation. There are many ways to do this, but we shall implement it by leveraging a maximum of BeepBeep functions, and minimizing the amount of custom code.

• We first define a BeepBeep Function object, whose task is to turn a Map produced by the ArenaSource into a Java BufferedImage. That is, a picture is created for each state of the arena; this is taken care of by the DrawArena class, which is made of a little more than 50 lines of code.
• These pictures must be shown somewhere; the BitmapJFrame is a very simple window object, containing a single JLabel widget whose background will be used to display the image. The window and its associated MouseListener make up 70 lines of code.

These 120 lines are the only pieces of custom code we need; the remainder of the program can be accomplished by creating and piping BeepBeep objects.

Among the various extensions (called "palettes") available for BeepBeep, one of them is called Widgets. In this palette, the WidgetSink processor has the task of receiving a stream of objects, and using them to update the state of a given Swing widget. For example, when the WidgetSink receives an image and is pointed to a JLabel, it will set the background property of the label to contain the image. Therefore, by simply piping the output of DrawArena to a WidgetSink that points to the label of our BitmapJFrame, the frame's content will automatically update every time a new image comes in.

Connecting the ArenaSource, DrawArena and WidgetSink together directly will produce nothing. This is caused by the fact that none of these processors produce events unless asked for. In order to drive the animation, a special BeepBeep processor, called the Pump, must be inserted into the chain. If inserted directly after the source, the pump will periodically query the arena for a new event (a process called "pulling"), and then "push" this new event into the drawing function and the remainder of the downstream chain. This final chain can be illustrated graphically:

In this drawing, events flow from left to right. It shows how the ArenaSource (leftmost box) is connected to a pump, itself connected to an ApplyFunction processor that applies the DrawArena function, which pipes its output into the WidgetSink responsible for displaying the image. The drawing also follows the BeepBeep graphical convention that each type of event is given a different color. Here, the maps produced by the arena are shown in pink, while the binary images produced by DrawArena are in light green. The Java code that creates this chain if made of just a few lines:

BitmapJFrame window = ...
Arena a = ...
JLabel l = window.getLabel();
ArenaSource as = new ArenaSource(a);
Pump pump = new Pump(50);
Connector.connect(as, pump);
ApplyFunction af = new ApplyFunction(new DrawArena(W, H));
Connector.connect(pump, af);
WidgetSink ws = new WidgetSink(l);
Connector.connect(af, ws);

Displaying the window and running the animation is just a matter of making the frame visible and starting the pump:

window.setVisible(true);
pump.start();

In our code sample, the pump was instantiated with parameter 50: this indicates that, once started, it will pull one new input event every 50 ms. This will result in an "animated" arena that updates at a rate of roughly 20 images per second. (In our program, the pump is actually started when the user clicks on the window.)

Plotting the evolution

The second part of the Post simulator is a dynamic plot that shows the evolution of the number of healthy, infected and recovered patients during the simulation. We shall again use BeepBeep's processors and functions to achieve a similar result, with special care to minimize the amount of custom code.

A single custom Function object needs to be created for this part of the program. The class GetHealth is a BeepBeep function that takes as input a Patient object, and returns as output the value of its Health status; it is made of 16 lines of code.

Equipped with this function, the remainder of the operation can be made using BeepBeep and its existing palettes. In order to do so, two palettes are needed:

• The Tuple palette provides basic functionalities for manipulating tuples (i.e. sets of key-value pairs)
• The MTNP palette allows users to create tables and display their contents by making background calls to Gnuplot

Let us start with the chain of processors, which we will then explain step by step.

The chain starts on the far left by a CountDecimate processor, which is instructed to keep only one out of every 25 input events. This is intended to reduce to lower the event rate in the remainder of the chain. The output is then divided into two copies. The top path turns every event into the number 1, and sends this feed of "ones" into a processor that adds them; this is a standard BeepBeep construct to create a counter (1, 2, 3, etc.) out of an arbitrary input stream. These values are then turned into a feed of tuples, by assigning them to an attribute named "t".

The bottom path takes the map of players in a given state, and applies a stack of functions to it, which reads from bottom to top. First, the map is transformed so that its values become the health status of each player; this is done by applying our custom GetHealth function to each value. The multiset of values in the map is then computed (function labelled **), the cardinality of each value is extracted (function labelled ##), and this map of value/cardinality pairs is finally turned into a tuple.

Both the top and the bottom paths produce a tuple; these two tuples are then merged using the "union" function. The end result of this first part of the chain is a stream of tuples, each made of t (the increasing timestamp counter) and one key-value pair each for the number of HEALTHY, INFECTED, RECOVERED and DEAD patients in the original map.

These tuples are then progressively accumulated into a table, that is sent to a DrawPlot processor producing a line plot out of its contents. The resulting image (produced by Gnuplot in the background) is then sent to a WidgetSink to be displayed in a window, as in the first part of our program.

The final result of this chain, which can then be connected to the ArenaSource, is a second window that displays the dynamically updated plot of the number of healthy, infected and recovered patients over time. Apart from the custom GetHealth function, creating this chain only requires piping existing BeepBeep objects and takes fewer than 25 Java instructions.

Wrapping up

We have used the current events about COVID-19 as a "fun" pretext to showcase some of the features of two important libraries developed at LIF.

• Complex data structures and objects can be synthesized according to various parameters using the Synthia generation library.
• The piping and processing of events can be handled by the BeepBeep library; in the present case, chains of a handful of basic processors and functions made it possible to display an animated simulation, as well as a dynamic plot computed from the sequence of states of that simulation.

In addition, we have seen how, using these libraries, the creation of the simulator can be done using a very limited amount of custom classes and lines of code.

This small project could lend itself to multiple easy extensions. Further use of Synthia's Picker objects could be used to generate more complex initial states or behaviors for each patient (e.g.: the possibility of being reinfected, or the separation of the infected state into asymptomatic and symptomatic, etc.). BeepBeep could also be used to perform further computations on the state of the simulation, such as the number of new infections over a sliding window (we leave this as an exercise!).

Running this program

To compile and run this example, first get the latest build of BeepBeep and its MTNP, Tuples, and Widgets palettes, and make sure they are in the classpath. The project also comes with a standardized Ant build script that has its own documentation.

In order to display the plots, Gnuplot must be installed and accessible by running gnuplot from the command line.

About the author

This example was coded by Sylvain Hallé, Full Professor at Université du Québec à Chicoutimi and Canada Research Chair on Software Testing, Specification and Verification.