0xDEADBEEF

covid-simulation.d

import std.algorithm;
import std.array;
import std.conv;
import std.getopt;
import std.random;
import std.range;
import std.stdio;

enum ushort Dim = 1000;

struct XY { int x; int y; }

struct Person {
  int infectionStep = -1;
  XY[] contacts;
  bool infected() { return infectionStep != -1; }
}

auto randomPos() {
  return XY(uniform(0, Dim), uniform(0, Dim));
}

XY[] generateContacts(XY xy, int contacts, int radius, int regionSize, double separationStrength) {
  return sequence!((_, i) => genContact(xy, radius, regionSize, separationStrength)).take(contacts).array;
}

XY genContact(XY xy, int radius, int regionSize, double separationStrength) {
  if (uniform01 < separationStrength) {
    XY regCorner1 = XY(xy.x - xy.x%regionSize, xy.y - xy.y%regionSize);
    XY regCorner2 = XY(regCorner1.x + regionSize, regCorner1.y + regionSize);
    return XY(
      uniform(regCorner1.x, regCorner2.x),
      uniform(regCorner1.y, regCorner2.y)
    );
  } else {
    return XY(
      uniform(max(0, xy.x-radius), min(Dim, xy.x+radius)),
      uniform(max(0, xy.y-radius), min(Dim, xy.y+radius))
    );
  }
}

// day 0 - infection
// day 3 - infectious
// day 5 - syptoms
bool isInfectious(ushort currentStep, ref Person p) {
  return p.infected && currentStep >= p.infectionStep+2 && currentStep <= p.infectionStep+5;
}

auto tryToInfect(ushort currentStep, ref Person p, double prob) {
  if (!p.infected && uniform01() < prob) {
    p.infectionStep = currentStep;
    return true;
  }
  return false;
}



void main(string[] args) {
  auto radius = Dim;
  auto infectionProbability = 0.1;
  auto contacts = 10;
  auto regionSize = 10;
  auto separationStrength = 0.5;

  auto help = getopt(
    args,
    std.getopt.config.passThrough,
    "radius|r",   "radius from which contacts are generated", &radius,
    "prob|p",     "infection probability", &infectionProbability,
    "contacts|c", "number of contacts", &contacts,
    "regionSize|s", "", &regionSize,
    "separationStrength|x", "", &separationStrength
  );

  if (help.helpWanted) {
    defaultGetoptPrinter("options:", help.options);
    return;
  }

  writeln("radius: ", radius);
  writeln("infectionProbability: ", infectionProbability);
  writeln("contacts: ", contacts);
  writeln("regionSize: ", regionSize);
  writeln("separationStrength: ", separationStrength);



  auto people = new Person[Dim][Dim];

  ref auto peopleXY(XY xy) { return people[xy.x][xy.y]; }


  foreach (ushort x; 0 .. Dim) {
    foreach (ushort y; 0 .. Dim) {
      people[x][y] = Person(-1, generateContacts(XY(x, y), contacts, radius, regionSize, separationStrength));
    }
  }

  foreach (_; 0 .. 100) {
    peopleXY(randomPos()).infectionStep = 0;
  }

  ulong nInfected = 0;
  foreach (ushort step; 0 .. 50) {

    foreach (ushort x; 0 .. Dim) {
      foreach (ushort y; 0 .. Dim) {
        auto p = people[x][y];
        if (!isInfectious(step, p)) continue;

        foreach (contact; p.contacts) {
          nInfected += tryToInfect(step, peopleXY(contact), infectionProbability);
        }
      }
    }

    auto infectedPercent = double(nInfected) / (Dim*Dim);// * 100;
    writefln!"step: %d, infected: %f"(step, infectedPercent);
  }
}