Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #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 = "30_"; // 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() {
- 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;
- }
Add Comment
Please, Sign In to add comment