Dependence of population evolution on age of giving birth

population_evolution.cpp

#include <fstream>
#include <iostream>
#include <random>
#include <set>
#include <string>
#include <tuple>

// Birth/death events
// If birth: time, +1, 0
// If death: time, -1, birthday
typedef std::tuple<int, int, int> event;

int initial_population = 10000;    // Number of people born at t=0
std::multiset<int> birthdays;      // Birthdays of currently alive people
std::multiset<event> events;       // Future births/deaths to be processed
int simulation_time = 2000;        // Number of years to run simulation for
int output_time = 10;              // Frequency to output population pyramid
int histogram_width = 5;           // Bin width for population histogram
std::string output_prefix = "";    // Prefix for output files
int birth_mean = 0;                // Mean year of giving birth
int birth_sigma = 2;               // Standard deviation year of giving birth
int birth_cutoff = 15;             // Lower bound for year of giving birth
int death_start = 50;              // Age at which people start dying
double death_start_probability = 0.005; // Probability of dying in start year
double death_probability_increment =
    0.005; // Increase of probability of dying every year after

int next_output_time = 0;   // Next moment to output population pyramid
int current_population = 0; // Current population

// Random number generators
std::random_device rd{};
std::mt19937 gen{rd()};
std::normal_distribution<double> parenthood_age;      // Age of having the child
std::uniform_real_distribution<double> uniform(0, 1); // rand()

void output() {
  std::ofstream out(output_prefix + std::to_string(next_output_time) + ".dat");
  std::cout << next_output_time << " " << birthdays.size() << std::endl;

  // Population pyramid starts with start_age and increments by histogram_width
  int start_age = 0, num_people = 0;
  for (auto sit = birthdays.rbegin();; ++sit) {
    if (sit == birthdays.rend()) {
      out << start_age << " " << num_people << "\n";
      out << start_age + histogram_width << " " << 0 << "\n";
      break;
    }

    int current_age = next_output_time - *sit;
    if (current_age >= start_age + histogram_width) {
      // If the age is outside the bin, output the bin and reset
      out << start_age << " " << num_people << "\n";
      start_age += histogram_width;
      num_people = 0;
      --sit;
    } else {
      // Otherwise update the count in the bin
      ++num_people;
    }
  }
}

// Process a birth/death
bool process() {
  int event_time, event_type, birthday;
  std::tie(event_time, event_type, birthday) = *(events.begin());

  // Output population pyramid if necessary
  if (event_time > next_output_time) {
    output();
    next_output_time += output_time;
    return true;
  }
  events.erase(events.begin());

  // If we have reached the end of the simulation, say goodbye
  if (event_time >= simulation_time)
    return false;

  if (event_type == -1) {
    // Death
    // Reduce population and remove birthday
    --current_population;
    birthdays.erase(birthdays.find(birthday));
  } else {
    // Birth
    // Increase population and insert birthday
    ++current_population;
    birthdays.insert(event_time);

    // Find birthday of child and add birth event
    int parenthood_age;
    do {
      parenthood_age = (int) ::parenthood_age(gen);
    } while(parenthood_age < birth_cutoff or parenthood_age >= death_start);
    events.insert({event_time + parenthood_age, 1, 0});

    // Find deathday and add death event
    int death_age = death_start;
    while(1) {
      double random_value = uniform(gen);
      if (random_value <=
          death_start_probability +
              (death_age - death_start) * death_probability_increment)
        break;
      ++death_age;
    }
    events.insert({event_time + death_age, -1, event_time});
  }
  return true;
}

void initialise() {
  output_prefix = std::to_string(birth_mean) + "_";
  parenthood_age = std::normal_distribution<double>(birth_mean, birth_sigma);
  for (int i = 0; i < initial_population; ++i) {
    events.insert({0, 1, 0});
  }
}

int main() { 
  std::cin >> birth_mean;
  initialise();
  while(process()) {};
  return 0;
}

Made for the Quora answer: How does delayed marriage help in population management?