# Be an Epidemiologist!

## Challenge

You must create a simple model of how disease spreads around a group of people.

## Rules and Requirements

The model must be a 1000 by 1000 2D array with each element being a different person.

The user must input three variables using argv: the probability of transmission (how likely someone is going to infect someone else), chance of mutation and how many periods the simulation should run for.

In the first period (t=0), four people should be chosen randomly and infected with the disease.

The way the disease behaves is governed by the following rules:

• The disease can only move vertically and horizontally, moving to the person next door.
• The infection last for 3 periods in every person. You may not factor in immunodeficiencies.
• After a person has been infected three times, they are immune and cannot be infected again.
• The disease is subject to mutation making the previously immune people vulnerable to this new mutated disease. The mutated disease has exactly the same traits and follows the same rules as the original disease.
• If a mutation occurs, the whole disease does not change, just that particular 'packet' upon transmission.
• Once a person has been been infected by one virus, they cannot be infected again until the current infection passes.
• If a person is infected, they are infectious from the start of their infection period to its end.
• There are no levels of immunity - a person is either immune or not.
• To stop memory overload, there is a maximum limit of 800 mutations.

At the end of the specified number of periods, you should output the results.

The results must be a 1000 x 1000 grid showing which people are infected and which people are not. This can be output as a text file, as an image file or graphical output (where #FFFFFF is a healthy person and #40FF00 is an infected person).

## Winning

The fastest code to run on my computer wins. Its time will be measured with the following piece of Python code:

import time, os
start = time.time()
os.system(command)
end = time.time()
print(end-start)


Note that when running this script, I will use the following defaults:

Probability of transmission = 1
Chance of mutation = 0.01
Number of periods = 1000

• You want to make a 10-gigabyte file? – Ypnypn Oct 2 '14 at 17:07
• By having a 4 GB limit, you have completely removed the option of saving the output in image file ... – Optimizer Oct 2 '14 at 18:16
• 1000x1000: That's more like it! – COTO Oct 2 '14 at 18:25
• Also say there are two people next to each other. The first contracts virus V, the second contracts virus V'. The contraction will both end at the same period. Can virus V infects the second person? (Or a more black-and-white question: is it possible that a person be infected immediately after he is healed, so he will end up with 6 consecutive period of infection?) – justhalf Oct 3 '14 at 2:04
• Another one, can two independent viruses mutate to the same virus? Say we have V in person A, and V again in person B. When they transmit the virus, can they both mutate to the same mutation V'? Or perhaps they in fact should mutate to the same virus strain? If they can mutate arbitrarily, what is the probability of two viruses mutating to the same virus strain? – justhalf Oct 3 '14 at 2:06

I was curious what this would look like so I made this quick and dirty hack in JavaScript: http://jsfiddle.net/andrewmaxwell/r8m54t9c/

// The probability that a healthy cell will be infected by an adjacent infected cell EACH FRAME.
var infectionProbability = 0.2

// The probability that the infection will mutate on transmission EACH FRAME.
var mutationProbability = 0.00001

// The maximum number of times a cell can be infected by the same infection.
var maxInfections = 3

// The number of frames a cell in infected before it becomes healthy again.
var infectionDuration = 3

// The width and heigh of the board
var size = 400

// The number of cells infected at the beginning.
var startingNum = 4

var imageData, // the visual representation of the board
cells, // array of cells
infectionCount // counter that is incremented whenever a mutation occurs

// Just some colors. The colors are re-used as the number of mutations increases.
var colors = [[255,0,0],[255,255,0],[0,255,0],[0,255,255],[0,0,255],[255,0,255],[128,0,0],[128,128,0],[0,128,0],[0,128,128],[0,0,128],[128,0,128],[255,128,128],[255,255,128],[128,255,128],[128,255,255],[128,128,255],[255,128,255]
]

// when a cell is infected, it isn't contagious until the next frame
function infect(person, infection){
person.infect = true
person.infectionCounts[infection] = (person.infectionCounts[infection] || 0) + 1
person.currentInfection = infection
}

// when a mutation occurs, it is given a number and the counter is incremented
function mutation(){
return infectionCount++
}

function reset(){

cells = []
infectionCount = 0
imageData = T.createImageData(size, size)

// initialize the cells, store them in a grid temporarily and an array for use in each frame
var grid = []
for (var i = 0; i < size; i++){
grid[i] = []
for (var j = 0; j < size; j++){
cells.push(grid[i][j] = {
infectionTime: 0, // how many frames until they are no longer infected, so 0 is healthy
infectionCounts: [], // this stores how many times the cell has been infected by each mutation
neighbors: [] // the neighboring cells
})
}
}

// store the neighbors of each cell, I just want to minimize the work done each frame
var neighborCoords = [[0,-1],[1,0],[0,1],[-1,0]]
for (var i = 0; i < size; i++){
for (var j = 0; j < size; j++){
for (var n = 0; n < neighborCoords.length; n++){
var row = i + neighborCoords[n]
var col = j + neighborCoords[n]
if (grid[row] && grid[row][col]){
grid[i][j].neighbors.push(grid[row][col])
}
}
}
}

// infect the initial cells
for (var i = 0; i < startingNum; i++){
infect(cells[Math.floor(cells.length * Math.random())], 0)
}
}

function loop(){
requestAnimationFrame(loop)

// for each cell marked as infected, set its infectionTime
for (var i = 0; i < cells.length; i++){
var p = cells[i]
if (p.infect){
p.infect = false
p.infectionTime = infectionDuration
}
}

for (var i = 0; i < cells.length; i++){
var p = cells[i]

// for each infected cell, decrement its timer
if (p.infectionTime){
p.infectionTime--

// for each neighbor that isn't infected, if the probability is right and the neighbor isn't immune to that infection, infect it
for (var n = 0; n < p.neighbors.length; n++){
var neighbor = p.neighbors[n]
if (!neighbor.infectionTime && Math.random() < infectionProbability){
var infection = Math.random() < mutationProbability ? mutation() : p.currentInfection
if (!neighbor.infectionCounts[infection] || neighbor.infectionCounts[infection] < maxInfections){
infect(neighbor, infection)
}
}
}

// colors! yay!
var color = colors[p.currentInfection % colors.length]
imageData.data[4 * i + 0] = color
imageData.data[4 * i + 1] = color
imageData.data[4 * i + 2] = color
} else {
imageData.data[4 * i + 0] = imageData.data[4 * i + 1] = imageData.data[4 * i + 2] = 0
}

imageData.data[4 * i + 3] = 255
}

T.putImageData(imageData, 0, 0)
}

// init canvas and go
C.width = C.height = size
T = C.getContext('2d')
reset()
loop()

• Setting infectionProbability to 1 made some of the sweetest patterns I've seen! – William Barbosa Oct 17 '14 at 17:37

# C++11, 6-8 minutes

My test run takes about 6-8 minutes in my Fedora 19, i5 machine. But due to the randomness of the mutation, it might as well be faster or take longer than that. I think the scoring criteria needs to be readdressed.

Prints the result as text at the end of completion, healthy person denoted by dot (.), infected person by asterisk (*), unless ANIMATE flag is set to true, in which case it will display different characters for people infected with different virus strain.

Here is a GIF for 10x10, 200 periods. ## Mutation behaviour

Each mutation will give new strain never seen before (so it is possible that one person infects the four neighboring people with 4 distinct strains), unless 800 strains have been generated, in which case no virus will go any further mutation.

The 8-minute result comes from the following number of infected people:

Period   0, Infected:      4
Period 100, Infected:  53743
Period 200, Infected: 134451
Period 300, Infected: 173369
Period 400, Infected: 228176
Period 500, Infected: 261473
Period 600, Infected: 276086
Period 700, Infected: 265774
Period 800, Infected: 236828
Period 900, Infected: 221275


while the 6-minute result comes from the following:

Period   0, Infected:      4
Period 100, Infected:  53627
Period 200, Infected: 129033
Period 300, Infected: 186127
Period 400, Infected: 213633
Period 500, Infected: 193702
Period 600, Infected: 173995
Period 700, Infected: 157966
Period 800, Infected: 138281
Period 900, Infected: 129381


## Person representation

Each person is represented in 205 bytes. Four bytes to store the virus type this person is contracting, one byte to store how long has this person been infected, and 200 bytes to store how many times he has contracted each strain of virus (2 bit each). Perhaps there are some additional byte-alignment done by C++, but the total size will be around 200MB. I have two grids to store the next step, so total it uses around 400MB.

I store the location of infected people in a queue, to cut the time required in the early periods (which is really useful up to periods < 400).

## Program technicalities

Every 100 steps this program will print the number of infected people, unless ANIMATE flag is set to true, in which case it will print the whole grid every 100ms.

This requires C++11 libraries (compile using -std=c++11 flag, or in Mac with clang++ -std=c++11 -stdlib=libc++ virus_spread.cpp -o virus_spread).

Run it without arguments for the default values, or with arguments like this:

./virus_spread 1 0.01 1000

#include <cstdio>
#include <cstring>
#include <random>
#include <cstdlib>
#include <utility>
#include <iostream>
#include <deque>
#include <cmath>
#include <functional>
#include <unistd.h>

typedef std::pair<int, int> pair;
typedef std::deque<pair> queue;

const bool ANIMATE = false;
const int MY_RAND_MAX = 999999;

std::default_random_engine generator(time(0));
std::uniform_int_distribution<int> distInt(0, MY_RAND_MAX);
auto randint = std::bind(distInt, generator);
std::uniform_real_distribution<double> distReal(0, 1);
auto randreal = std::bind(distReal, generator);

const int VIRUS_TYPE_COUNT = 800;
const int SIZE = 1000;
const int VIRUS_START_COUNT = 4;

typedef struct Person{
int virusType;
char time;
uint32_t immune[VIRUS_TYPE_COUNT/16];
} Person;

Person people[SIZE][SIZE];
Person tmp[SIZE][SIZE];
queue infecteds;

double transmissionProb = 1.0;
double mutationProb = 0.01;
int periods = 1000;

char inline getTime(Person person){
return person.time;
}

char inline getTime(int row, int col){
return getTime(people[row][col]);
}

Person inline setTime(Person person, char time){
person.time = time;
return person;
}

Person inline addImmune(Person person, uint32_t type){
person.immune[type/16] += 1 << (2*(type % 16));
return person;
}

bool inline infected(Person person){
return getTime(person) > 0;
}

bool inline infected(int row, int col){
return infected(tmp[row][col]);
}

bool inline immune(Person person, uint32_t type){
return (person.immune[type/16] >> (2*(type % 16)) & 3) == 3;
}

bool inline immune(int row, int col, uint32_t type){
return immune(people[row][col], type);
}

Person inline infect(Person person, uint32_t type){
person.time = 1;
person.virusType = type;
return person;
}

bool inline infect(int row, int col, uint32_t type){
auto person = people[row][col];
auto tmpPerson = tmp[row][col];
if(infected(tmpPerson) || immune(tmpPerson, type) || infected(person) || immune(person, type)) return false;
person = infect(person, type);
infecteds.push_back(std::make_pair(row, col));
tmp[row][col] = person;
return true;
}

uint32_t inline getType(Person person){
return person.virusType;
}

uint32_t inline getType(int row, int col){
return getType(people[row][col]);
}

void print(){
for(int row=0; row < SIZE; row++){
for(int col=0; col < SIZE; col++){
printf("%c", infected(row, col) ? (ANIMATE ? getType(row, col)+48 : '*') : '.');
}
printf("\n");
}
}

void move(){
for(int row=0; row<SIZE; ++row){
for(int col=0; col<SIZE; ++col){
people[row][col] = tmp[row][col];
}
}
}

int main(const int argc, const char **argv){
if(argc > 3){
transmissionProb = std::stod(argv);
mutationProb = std::stod(argv);
periods = atoi(argv);
}
int row, col, size;
uint32_t type, newType=0;
char time;
Person person;
memset(people, 0, sizeof(people));
for(int row=0; row<SIZE; ++row){
for(int col=0; col<SIZE; ++col){
people[row][col] = {};
}
}
for(int i=0; i<VIRUS_START_COUNT; i++){
row = randint() % SIZE;
col = randint() % SIZE;
if(!infected(row, col)){
infect(row, col, 0);
} else {
i--;
}
}
move();
if(ANIMATE){
print();
}
for(int period=0; period < periods; ++period){
size = infecteds.size();
for(int i=0; i<size; ++i){
pair it = infecteds.front();
infecteds.pop_front();
row = it.first;
col = it.second;
person = people[row][col];
time = getTime(person);
if(time == 0) continue;
type = getType(person);
if(row > 0 && randreal() < transmissionProb){
if(newType < VIRUS_TYPE_COUNT-1 && randreal() < mutationProb){
newType++;
if(!infect(row-1, col, newType)) newType--;
} else {
infect(row-1, col, type);
}
}
if(row < SIZE-1 && randreal() < transmissionProb){
if(newType < VIRUS_TYPE_COUNT-1 && randreal() < mutationProb){
newType++;
if(!infect(row+1, col, newType)) newType--;
} else {
infect(row+1, col, type);
}
}
if(col > 0 && randreal() < transmissionProb){
if(newType < VIRUS_TYPE_COUNT-1 && randreal() < mutationProb){
newType++;
if(!infect(row, col-1, newType)) newType--;
} else {
infect(row, col-1, type);
}
}
if(col < SIZE-1 && randreal() < transmissionProb){
if(newType < VIRUS_TYPE_COUNT-1 && randreal() < mutationProb){
newType++;
if(!infect(row, col+1, newType)) newType--;
} else {
infect(row, col+1, type);
}
}
time += 1;
if(time == 4) time = 0;
person = setTime(person, time);
if(time == 0){
} else {
infecteds.push_back(std::make_pair(row, col));
}
tmp[row][col] = person;
}
if(!ANIMATE && period % 100 == 0) printf("Period %d, Size: %d\n", period, size);
move();
if(ANIMATE){
printf("\n");
print();
usleep(100000);
}
}
if(!ANIMATE){
print();
}
return 0;
}

• I really like this! My only question is how do you make the GIF? – Beta Decay Oct 4 '14 at 8:56
• I'm using this tool: linux.die.net/man/1/byzanz-record . Currently it doesn't have GUI, so you'll need to use the command line =D – justhalf Oct 4 '14 at 9:56
• Oh, that's nice, thanks! :) – Beta Decay Oct 4 '14 at 10:04

# Edit 2

I've finally(5 hours) generated a verbose output for close to 1000 periods(only 840 frames then it crashed) at 1000x1000, every 1 period, however its close to 160MB and requires all the memory on my system to display(IrfanView), not even sure that would work in a browser, may I'll put it up later.

# EDIT

I've spent allot of time on making this more efficient per "Beta Decay"s answer stating "Choose the strain randomly" I have choosin only the random method for choosing who infects who per period, however I've changed the way that's calculated and threaded out every thing, I've updated my posts with the new details.

Coded my closest estimation to this I could, I hope it follows all the rules, it uses a ton of memory on my system(about 1.2GB). The program can output animated gifs(looks cool, really slow) or just a image matching "Beta Decay"s specs. This is a little bit of reinventing the wheel, but definitely looks cool:

# Results

(Note: this only differentiates between infected and non infected, ie non-verbose)

## 1000 Periods, 1% Mutation Rate, 100% Spread: # Examples(Verbose)

Anyway using 100% "Probability of transmission" in non-verbose mode is kind of boring as you always get the same shapes and you can't see the different mutations, if you tweak the parameters around a-bit(and enable verbose mode) you get some cool looking output(animated GIFs display every 10th frame):

## Random - Grid Size: 200, ProbTransmission: 100%, ProbMutation: 1% ## Random - Grid Size: 200, ProbTransmission: 20%, ProbMutation: 1% # Scoring

I agree with "justhalf" that the scoring criteria may not be fair as every run will differ due to the randomness of mutations and position of random start points. Maybe we could do average of several runs or something like that..., well anyway this is in C# so bounty for me :( anyway.

# Code

Make sure to include the MagickImage library(set to compile x64 bit) otherwise it won't build(http://pastebin.com/vEmPF1PM):

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
using System.Drawing.Imaging;
using ImageMagick;
using System.IO;

namespace Infection
{
class Program
{
#region Infection Options
private const double ProbabilityOfTransmission = .2;
private const double ChanceOfMutation = 0.01;
private const Int16 StageSize = 1000;
private const Int16 MaxNumberOfMutations = 800;
private const byte MaxInfectionTime = 3;
private const byte NumberOfPeopleToRandomlyInfect = 4;
private static int NumberOfPeriods = 1000;
#endregion Infection Options

#region Run Options
private const bool VerbosMode = false;
private const int ImageFrequency = 10;
#endregion Run Options

#region Stage
private static Int16 MutationNumber = 1;

private class Person
{
public Person()
{
PreviousInfections = new Dictionary<Int16, byte>();
InfectionTime = 0;
CurrentInfection = 0;
PossibleNewInfections = new List<short>(4);
}
public Dictionary<Int16, byte> PreviousInfections { get; set; }
public byte InfectionTime { get; set; }
public Int16 CurrentInfection { get; set; }
public List<Int16> PossibleNewInfections { get; set; }
}
private static Person[][] Stage = new Person[StageSize][];
#endregion Stage

static void Main(string[] args)
{
DateTime start = DateTime.UtcNow;

//Initialize stage
for (Int16 i = 0; i < Stage.Length; i++)
{
var tmpList = new List<Person>();
for (Int16 j = 0; j < Stage.Length; j++)
Stage[i] = tmpList.ToArray();
}

//Randomly infect people
RandomlyInfectPeople(NumberOfPeopleToRandomlyInfect);

//Run through the periods(NumberOfPeriods times)
List<MagickImage> output = new List<MagickImage>();
while (NumberOfPeriods > 0)
{
//Print details(verbose)
if (VerbosMode && NumberOfPeriods % ImageFrequency == 0)
{
Console.WriteLine("Current Number: " + NumberOfPeriods);
Console.WriteLine("Current Mutation: " + MutationNumber);
}

Period();
}

//Outputs a Animated Gif(verbose)
if (VerbosMode)
{
ImagesToAnimatedGIF(output.ToArray(), Directory.GetCurrentDirectory() + "\\Output.gif");
System.Diagnostics.Process.Start(Directory.GetCurrentDirectory() + "\\Output.gif");
}
//Only outputs the basic result image matching the specs
SaveBoardToSimpleImage(Directory.GetCurrentDirectory() + "\\FinalState.gif");

Console.WriteLine("Total run time in seconds: " + (DateTime.UtcNow - start).TotalSeconds);
Console.WriteLine("Press enter to exit");
}

#region Image
private static void SaveBoardToSimpleImage(string filepath)
{
using (Bitmap img = new Bitmap(StageSize, StageSize))
{
for (int i = 0; i < img.Width; i++)
for (int j = 0; j < img.Height; j++)
img.SetPixel(i, j, Stage[i][j].CurrentInfection == 0 ? Color.FromArgb(255, 255, 255) :
Color.FromArgb(64, 255, 0));
img.Save(filepath, ImageFormat.Gif);
}
}
private static MagickImage BoardToImage()
{
using (Bitmap img = new Bitmap(StageSize, StageSize))
{
for (int i = 0; i < img.Width; i++)
for (int j = 0; j < img.Height; j++)
img.SetPixel(i, j, Stage[i][j].CurrentInfection == 0 ? Color.White :
Color.FromArgb(Stage[i][j].CurrentInfection % 255,
Math.Abs(Stage[i][j].CurrentInfection - 255) % 255,
Math.Abs(Stage[i][j].CurrentInfection - 510) % 255));
return new MagickImage(img);
}
}
private static void ImagesToAnimatedGIF(MagickImage[] images, string filepath)
{
using (MagickImageCollection collection = new MagickImageCollection())
{
foreach (var image in images)
{
collection.Last().AnimationDelay = 20;
}
collection.Write(filepath);
}
}
#endregion Image

#region Infection
private static void Period()
{
Infect();
ChooseRandomInfections();
IncrementDiseaseProgress();
Cure();

NumberOfPeriods--;
}
private static void Cure()
{
Parallel.For(0, Stage.Length, i =>
{
for (Int16 j = 0; j < Stage.Length; j++)
if (Stage[i][j].CurrentInfection != 0 && Stage[i][j].InfectionTime == MaxInfectionTime + 1)
{
if (Stage[i][j].PreviousInfections.ContainsKey(Stage[i][j].CurrentInfection))
Stage[i][j].PreviousInfections[Stage[i][j].CurrentInfection]++;
else

//Cure
Stage[i][j].InfectionTime = 0;
Stage[i][j].CurrentInfection = 0;
}
});
}
private static void IncrementDiseaseProgress()
{
Parallel.For(0, Stage.Length, i =>
{
for (Int16 j = 0; j < Stage.Length; j++)
if (Stage[i][j].CurrentInfection != 0)
Stage[i][j].InfectionTime++;
});
}
private static void RandomlyInfectPeople(Int16 numberOfPeopleToInfect)
{
var randomList = new List<int>();
while (randomList.Count() < numberOfPeopleToInfect * 2)
{
randomList = randomList.Distinct().ToList();
}
while (randomList.Count() > 0)
{
Stage[randomList.Last()][randomList[randomList.Count() - 2]].CurrentInfection = MutationNumber;
Stage[randomList.Last()][randomList[randomList.Count() - 2]].InfectionTime = 1;
randomList.RemoveAt(randomList.Count() - 2);
randomList.RemoveAt(randomList.Count() - 1);
}
}
private static void Infect()
{
Parallel.For(0, Stage.Length, i =>
{
for (Int16 j = 0; j < Stage.Length; j++)
InfectAllSpacesAround((short)i, j);
});
}
private static void InfectAllSpacesAround(Int16 x, Int16 y)
{
//If not infected or just infected this turn return
if (Stage[x][y].CurrentInfection == 0 || (Stage[x][y].CurrentInfection != 0 && Stage[x][y].InfectionTime == 0)) return;

//Infect all four directions(if possible)
if (x > 0)
InfectOneSpace(Stage[x][y].CurrentInfection, (short)(x - 1), y);

if (x < Stage.Length - 1)
InfectOneSpace(Stage[x][y].CurrentInfection, (short)(x + 1), y);

if (y > 0)
InfectOneSpace(Stage[x][y].CurrentInfection, x, (short)(y - 1));

if (y < Stage.Length - 1)
InfectOneSpace(Stage[x][y].CurrentInfection, x, (short)(y + 1));
}
private static void InfectOneSpace(Int16 currentInfection, Int16 x, Int16 y)
{
//If the person is infected, or If they've already been infected "MaxInfectionTime" then don't infect
if (Stage[x][y].CurrentInfection != 0 || (Stage[x][y].PreviousInfections.ContainsKey(currentInfection) &&
Stage[x][y].PreviousInfections[currentInfection] >= MaxInfectionTime)) return;

//If random is larger than change of transmission don't transmite disease
if (RandomGen2.Next(100) + 1 > ProbabilityOfTransmission * 100) return;

//Possible mutate
if (MutationNumber <= MaxNumberOfMutations && RandomGen2.Next(100) + 1 <= ChanceOfMutation * 100)
lock (Stage[x][y])
{
MutationNumber++;
}
//Regular infection
else
lock (Stage[x][y])

}
private static void ChooseRandomInfections()
{
Parallel.For(0, Stage.Length, i =>
{
for (Int16 j = 0; j < Stage.Length; j++)
{
if (Stage[i][j].CurrentInfection != 0 || !Stage[i][j].PossibleNewInfections.Any()) continue;
Stage[i][j].CurrentInfection = Stage[i][j].PossibleNewInfections[RandomGen2.Next(Stage[i][j].PossibleNewInfections.Count)];
Stage[i][j].PossibleNewInfections.Clear();
Stage[i][j].InfectionTime = 0;
}
}
);
}
#endregion Infection
}

//Fancy Schmancy new random number generator for threaded stuff, fun times
//http://blogs.msdn.com/b/pfxteam/archive/2009/02/19/9434171.aspx
public static class RandomGen2
{
private static Random _global = new Random();
private static Random _local;

public static int Next()
{
Random inst = _local;
if (inst == null)
{
int seed;
lock (_global) seed = _global.Next();
_local = inst = new Random(seed);
}
return inst.Next();
}

public static int Next(int input)
{
Random inst = _local;
if (inst == null)
{
int seed;
lock (_global) seed = _global.Next();
_local = inst = new Random(seed);
}
return inst.Next(input);
}
}
}