Guest User

Untitled

a guest
Nov 17th, 2018
114
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 4.01 KB | None | 0 0
  1. #include <fstream>
  2. #include <iostream>
  3. #include <random>
  4. #include <set>
  5. #include <string>
  6. #include <tuple>
  7.  
  8. // Birth/death events
  9. // If birth: time, +1, 0
  10. // If death: time, -1, birthday
  11. typedef std::tuple<int, int, int> event;
  12.  
  13. int initial_population = 10000; // Number of people born at t=0
  14. std::multiset<int> birthdays; // Birthdays of currently alive people
  15. std::multiset<event> events; // Future births/deaths to be processed
  16. int simulation_time = 2000; // Number of years to run simulation for
  17. int output_time = 10; // Frequency to output population pyramid
  18. int histogram_width = 5; // Bin width for population histogram
  19. std::string output_prefix = "30_"; // Prefix for output files
  20. int birth_mean = 0; // Mean year of giving birth
  21. int birth_sigma = 2; // Standard deviation year of giving birth
  22. int birth_cutoff = 15; // Lower bound for year of giving birth
  23. int death_start = 50; // Age at which people start dying
  24. double death_start_probability = 0.005; // Probability of dying in start year
  25. double death_probability_increment =
  26. 0.005; // Increase of probability of dying every year after
  27.  
  28. int next_output_time = 0; // Next moment to output population pyramid
  29. int current_population = 0; // Current population
  30.  
  31. // Random number generators
  32. std::random_device rd{};
  33. std::mt19937 gen{rd()};
  34. std::normal_distribution<double> parenthood_age; // Age of having the child
  35. std::uniform_real_distribution<double> uniform(0, 1); // rand()
  36.  
  37. void output() {
  38. std::ofstream out(output_prefix + std::to_string(next_output_time) + ".dat");
  39. std::cout << next_output_time << " " << birthdays.size() << std::endl;
  40.  
  41. // Population pyramid starts with start_age and increments by histogram_width
  42. int start_age = 0, num_people = 0;
  43. for (auto sit = birthdays.rbegin();; ++sit) {
  44. if (sit == birthdays.rend()) {
  45. out << start_age << " " << num_people << "\n";
  46. out << start_age + histogram_width << " " << 0 << "\n";
  47. break;
  48. }
  49.  
  50. int current_age = next_output_time - *sit;
  51. if (current_age >= start_age + histogram_width) {
  52. // If the age is outside the bin, output the bin and reset
  53. out << start_age << " " << num_people << "\n";
  54. start_age += histogram_width;
  55. num_people = 0;
  56. --sit;
  57. } else {
  58. // Otherwise update the count in the bin
  59. ++num_people;
  60. }
  61. }
  62. }
  63.  
  64. // Process a birth/death
  65. bool process() {
  66. int event_time, event_type, birthday;
  67. std::tie(event_time, event_type, birthday) = *(events.begin());
  68.  
  69. // Output population pyramid if necessary
  70. if (event_time > next_output_time) {
  71. output();
  72. next_output_time += output_time;
  73. return true;
  74. }
  75. events.erase(events.begin());
  76.  
  77. // If we have reached the end of the simulation, say goodbye
  78. if (event_time >= simulation_time)
  79. return false;
  80.  
  81. if (event_type == -1) {
  82. // Death
  83. // Reduce population and remove birthday
  84. --current_population;
  85. birthdays.erase(birthdays.find(birthday));
  86. } else {
  87. // Birth
  88. // Increase population and insert birthday
  89. ++current_population;
  90. birthdays.insert(event_time);
  91.  
  92. // Find birthday of child and add birth event
  93. int parenthood_age;
  94. do {
  95. parenthood_age = (int) ::parenthood_age(gen);
  96. } while(parenthood_age < birth_cutoff or parenthood_age >= death_start);
  97. events.insert({event_time + parenthood_age, 1, 0});
  98.  
  99. // Find deathday and add death event
  100. int death_age = death_start;
  101. while(1) {
  102. double random_value = uniform(gen);
  103. if (random_value <=
  104. death_start_probability +
  105. (death_age - death_start) * death_probability_increment)
  106. break;
  107. ++death_age;
  108. }
  109. events.insert({event_time + death_age, -1, event_time});
  110. }
  111. return true;
  112. }
  113.  
  114. void initialise() {
  115. parenthood_age = std::normal_distribution<double>(birth_mean, birth_sigma);
  116. for (int i = 0; i < initial_population; ++i) {
  117. events.insert({0, 1, 0});
  118. }
  119. }
  120.  
  121. int main() {
  122. std::cin >> birth_mean;
  123. initialise();
  124. while(process()) {};
  125. return 0;
  126. }
Add Comment
Please, Sign In to add comment