Stride Reference Manual  1.0
PopulationGenerator.cpp
Go to the documentation of this file.
1 #include "PopulationGenerator.h"
2 #include "FamilyParser.h"
3 #include "util/InstallDirs.h"
4 #include "util/TimeStamp.h"
5 
6 #include <stdexcept>
7 #include <boost/property_tree/xml_parser.hpp>
8 #include <limits>
9 #include <algorithm>
10 
11 using namespace stride;
12 using namespace popgen;
13 using namespace util;
14 using namespace boost::property_tree;
15 using namespace xml_parser;
16 
17 template<class U>
18 PopulationGenerator<U>::PopulationGenerator(const string& filename, const int& seed, bool output) {
19  // check data environment.
20  if (InstallDirs::getDataDir().empty()) {
21  throw runtime_error(string(__func__) + "> Data directory not present! Aborting.");
22  }
23 
24  try {
25  read_xml((InstallDirs::getDataDir() /= filename).string(), m_props, trim_whitespace | no_comments);
26  } catch (exception& e) {
27  throw invalid_argument(string("Invalid file: ") + (InstallDirs::getDataDir() /= filename).string());
28  }
29 
30  try {
31  long int tot = m_props.get < long
32  int > ("population.<xmlattr>.total");
33 
34  if (tot <= 0) {
35  throw invalid_argument("Invalid attribute population::total");
36  }
37 
38  m_total = tot;
39  } catch (invalid_argument& e) {
40  throw e;
41  } catch (exception& e) {
42  throw invalid_argument("Missing/invalid attribute population::total");
43  }
44 
45  m_next_id = 1;
46  m_output = output;
47  m_rng = U(seed);
48 
49  checkForValidXML();
50 }
51 
52 template<class U>
53 void PopulationGenerator<U>::generate(const string& prefix) {
54 
55  if (m_output) cerr << "Generating " << m_total << " people...\n";
56 
57  makeHouseholds();
58  makeCities();
59  makeVillages();
60  placeHouseholds();
61  makeSchools();
62  makeUniversities();
63  makeWork();
64  makeCommunities();
65  assignToSchools();
66  assignToUniversities();
67  assignToWork();
68 
69  double start_radius = m_props.get<double>("population.commutingdata.<xmlattr>.start_radius");
70  double factor = m_props.get<double>("population.commutingdata.<xmlattr>.factor");
71 
72  vector<pair<GeoCoordinate, map<double, vector<uint>>>> distance_map = makeDistanceMap(start_radius, factor,
73  m_primary_communities);
74  assignToCommunities(distance_map, m_primary_communities, &SimplePerson::m_primary_community, "primary communities");
75 
76  distance_map.clear();
77  distance_map = makeDistanceMap(start_radius, factor, m_secondary_communities);
78  assignToCommunities(distance_map, m_secondary_communities, &SimplePerson::m_secondary_community,
79  "secondary communities");
80 
81  if (m_output) cerr << "Generated " << m_people.size() << " people\n";
82 
83  string target_cities = prefix + "_districts.csv";
84  string target_pop = prefix + "_people.csv";
85  string target_households = prefix + "_households.csv";
86  string target_clusters = prefix + "_clusters.csv";
87  string target_summary = prefix + ".xml";
88 
89  writeCities(target_cities);
90  writePop(target_pop);
91  writeHouseholds(target_households);
92  writeClusters(target_clusters);
93 
94  // Now write a summary
95  ptree config;
96  config.put("population.people", target_pop);
97  config.put("population.districts", target_cities);
98  config.put("population.clusters", target_clusters);
99  config.put("population.households", target_households);
100  config.add_child("population.cities", m_props.get_child("population.cities"));
101  write_xml((InstallDirs::getDataDir() /= target_summary).string(), config);
102  if (m_output) cout << "Written summary " << target_summary << endl;
103 }
104 
105 template<class U>
106 void PopulationGenerator<U>::writeCities(const string& target_cities) {
107  ofstream my_file {(InstallDirs::getDataDir() /= target_cities).string()};
108  double total_pop = 0.0;
109 
110  for (const SimpleCity& city: m_cities) {
111  total_pop += city.m_current_size;
112  }
113 
114  for (const SimpleCluster& village: m_villages) {
115  total_pop += village.m_current_size;
116  }
117 
118  if (my_file.is_open()) {
119  my_file
120  << "\"city_id\",\"city_name\",\"province\",\"population\",\"x_coord\",\"y_coord\",\"latitude\",\"longitude\"\n";
121 
122  uint provinces = m_props.get<uint>("population.<xmlattr>.provinces");
123  AliasDistribution dist {vector<double>(provinces, 1.0 / provinces)};
124 
125  auto printCityData = [&](const SimpleCity& to_print) {
126  my_file << to_print.m_current_size / total_pop
127  << ",0,0,"
128  << to_print.m_coord.m_latitude
129  << ","
130  << to_print.m_coord.m_longitude
131  << endl;
132  };
133 
134  auto printVillageData = [&](const SimpleCluster& to_print) {
135  SimpleCity city = SimpleCity(to_print.m_current_size, to_print.m_max_size, to_print.m_id, "",
136  to_print.m_coord);
137  printCityData(city);
138  };
139 
140  for (const SimpleCity& city: m_cities) {
141  my_file.precision(std::numeric_limits<double>::max_digits10);
142  my_file << city.m_id
143  << ",\""
144  << city.m_name
145  << "\"," << dist(m_rng) + 1 << ",";
146  printCityData(city);
147  }
148 
149  uint village_counter = 1;
150  for (const SimpleCluster& village: m_villages) {
151  my_file.precision(std::numeric_limits<double>::max_digits10);
152  my_file << village.m_id
153  << ",\""
154  << village_counter
155  << "\"," << dist(m_rng) + 1 << ",";
156 
157  printVillageData(village);
158  village_counter++;
159  }
160 
161  my_file.close();
162  if (m_output) cout << "Written " << target_cities << endl;
163  } else {
164  throw invalid_argument("In PopulationGenerator: Invalid file.");
165  }
166 }
167 
168 template<class U>
169 void PopulationGenerator<U>::writePop(const string& target_pop) const {
170  ofstream my_file {(InstallDirs::getDataDir() /= target_pop).string()};
171 
172  if (my_file.is_open()) {
173  my_file << "\"age\",\"household_id\",\"school_id\",\"work_id\",\"primary_community\",\"secondary_community\"\n";
174 
175  for (const SimplePerson& person: m_people) {
176  my_file << person.m_age << ","
177  << person.m_household_id << ","
178  << person.m_school_id << ","
179  << person.m_work_id << ","
180  << person.m_primary_community << ","
181  << person.m_secondary_community
182  << endl;
183  }
184 
185  my_file.close();
186  if (m_output) cout << "Written " << target_pop << endl;
187  } else {
188  throw invalid_argument("In PopulationGenerator: Invalid file.");
189  }
190 }
191 
192 template<class U>
193 void PopulationGenerator<U>::writeHouseholds(const string& target_households) const {
194  ofstream my_file {(InstallDirs::getDataDir() /= target_households).string()};
195  if (my_file.is_open()) {
196  my_file << "\"hh_id\",\"latitude\",\"longitude\",\"size\"\n";
197 
198  for (const SimpleHousehold& household: m_households) {
199  my_file << household.m_id << ","
200  << m_people.at(household.m_indices.at(0)).m_coord.m_latitude << ","
201  << m_people.at(household.m_indices.at(0)).m_coord.m_longitude << ","
202  << household.m_indices.size()
203  << endl;
204  }
205 
206  my_file.close();
207  if (m_output) cout << "Written " << target_households << endl;
208  } else {
209  throw invalid_argument("In PopulationGenerator: Invalid file.");
210  }
211 }
212 
213 template<class U>
214 void PopulationGenerator<U>::writeClusters(const string& target_clusters) const {
215  ofstream my_file {(InstallDirs::getDataDir() /= target_clusters).string()};
216  if (my_file.is_open()) {
217  my_file << "\"cluster_id\",\"cluster_type\",\"latitude\",\"longitude\"\n";
218 
219  vector<ClusterType> types {ClusterType::Household,
225 
226  for (auto& cluster_type: types) {
227  uint current_id = 1;
228  while (m_locations.find(make_pair(cluster_type, current_id)) != m_locations.end()) {
229  my_file.precision(std::numeric_limits<double>::max_digits10);
230  my_file << current_id << ","
231  << toString(cluster_type) << ","
232  << m_locations.at(make_pair(cluster_type, current_id)).m_latitude << ","
233  << m_locations.at(make_pair(cluster_type, current_id)).m_longitude
234  << endl;
235 
236  ++current_id;
237  }
238  }
239 
240  my_file.close();
241  if (m_output) cout << "Written " << target_clusters << endl;
242  } else {
243  throw invalid_argument("In PopulationGenerator: Invalid file.");
244  }
245 }
246 
247 template<class U>
249  ptree pop_config = m_props.get_child("population");
250 
253 
255  int provinces = pop_config.get<int>("<xmlattr>.provinces");
256 
257  if (provinces <= 0) {
258  throw invalid_argument("In PopulationGenerator: Numerical error.");
259  }
260 
262  double radius = pop_config.get<double>("commutingdata.<xmlattr>.start_radius");
263  double factor = pop_config.get<double>("commutingdata.<xmlattr>.factor");
264 
265  if (radius <= 0 || factor <= 1.0) {
266  throw invalid_argument("In PopulationGenerator: Numerical error.");
267  }
268 
270  if (m_output) cerr << "\rChecking for valid XML [0%]";
271  int total_size = 0;
272  bool has_no_cities = true;
273  auto cities_config = pop_config.get_child("cities");
274  vector<GeoCoordinate> current_locations;
275  for (auto it = cities_config.begin(); it != cities_config.end(); it++) {
276  if (it->first == "city") {
277  has_no_cities = false;
278  it->second.get<string>("<xmlattr>.name");
279  total_size += it->second.get<int>("<xmlattr>.pop");
280  double latitude = it->second.get<double>("<xmlattr>.lat");
281  double longitude = it->second.get<double>("<xmlattr>.lon");
282 
283  if (abs(latitude) > 90 || abs(longitude) > 180) {
284  throw invalid_argument("In PopulationGenerator: Invalid geo-coordinate in XML.");
285  }
286 
287  if (it->second.get<int>("<xmlattr>.pop") <= 0) {
288  throw invalid_argument("In PopulationGenerator: Numerical error.");
289  }
290 
291  auto it2 = find(current_locations.begin(), current_locations.end(), GeoCoordinate(latitude, longitude));
292  if (it2 != current_locations.end())
293  throw invalid_argument("In PopulationGenerator: Duplicate coordinates given in XML.");
294 
295  current_locations.push_back(GeoCoordinate(latitude, longitude));
296  } else {
297  throw invalid_argument("In PopulationGenerator: Missing/incorrect tags/attributes in XML.");
298  }
299  }
300 
301  if (has_no_cities) {
302  throw invalid_argument("In PopulationGenerator: No cities found.");
303  }
304 
306  if (m_output) cerr << "\rChecking for valid XML [18%]";
307  auto village_config = pop_config.get_child("villages");
308  double village_radius_factor = village_config.get<double>("<xmlattr>.radius");
309 
310  bool has_no_villages = true;
311 
312  double fraction = 0.0;
313  for (auto it = village_config.begin(); it != village_config.end(); it++) {
314  if (it->first == "village") {
315  has_no_villages = false;
316  int minimum = it->second.get<int>("<xmlattr>.min");
317  int max = it->second.get<int>("<xmlattr>.max");
318  fraction += it->second.get<double>("<xmlattr>.fraction") / 100.0;
319 
320  if (fraction < 0) {
321  throw invalid_argument("In PopulationGenerator: Numerical error.");
322  }
323 
324  if (minimum > max || minimum <= 0 || max < 0) {
325  throw invalid_argument("In PopulationGenerator: Numerical error.");
326  }
327  } else if (it->first == "<xmlattr>") {
328  } else {
329  throw invalid_argument("In PopulationGenerator: Missing/incorrect tags/attributes in XML.");
330  }
331  }
332 
333  if (fraction != 1.0 || village_radius_factor <= 0.0) {
334  throw invalid_argument("In PopulationGenerator: Numerical error.");
335  }
336 
337  if (has_no_villages) {
338  throw invalid_argument("In PopulationGenerator: No villages found.");
339  }
340 
341 
344  if (m_output) cerr << "\rChecking for valid XML [36%]";
345  auto education_config = pop_config.get_child("education");
346  auto school_work_config = pop_config.get_child("school_work_profile");
347  total_size = education_config.get<int>("mandatory.<xmlattr>.total_size");
348  int cluster_size = education_config.get<int>("mandatory.<xmlattr>.cluster_size");
349  int mandatory_min = school_work_config.get<int>("mandatory.<xmlattr>.min");
350  int mandatory_max = school_work_config.get<int>("mandatory.<xmlattr>.max");
351 
352  if (mandatory_min > mandatory_max || mandatory_min < 0 || mandatory_max < 0) {
353  throw invalid_argument("In PopulationGenerator: Numerical error in min/max pair.");
354  }
355 
356  if (total_size <= 0 || cluster_size <= 0) {
357  throw invalid_argument("In PopulationGenerator: Numerical error.");
358  }
359 
361  if (m_output) cerr << "\rChecking for valid XML [42%]";
362  school_work_config = pop_config.get_child("school_work_profile.employable.young_employee");
363  int minimum = school_work_config.get<int>("<xmlattr>.min");
364  int max = school_work_config.get<int>("<xmlattr>.max");
365  cluster_size = education_config.get<int>("optional.<xmlattr>.cluster_size");
366  fraction = 1.0 - school_work_config.get<double>("<xmlattr>.fraction") / 100.0;
367  total_size = education_config.get<uint>("optional.<xmlattr>.total_size");
368 
369  if (minimum > max || minimum < 0 || max < 0) {
370  throw invalid_argument("In PopulationGenerator: Numerical error in min/max pair.");
371  }
372 
373  if (total_size <= 0 || cluster_size <= 0 || fraction < 0.0 || fraction > 1.0) {
374  throw invalid_argument("In PopulationGenerator: Numerical error.");
375  }
376 
377  if (minimum <= mandatory_max) {
378  throw invalid_argument("In PopulationGenerator: Overlapping min/max pairs.");
379  }
380 
381  fraction = education_config.get<double>("optional.far.<xmlattr>.fraction") / 100.0;
382 
383  if (fraction < 0.0 || fraction > 1.0) {
384  throw invalid_argument("In PopulationGenerator: Numerical error.");
385  }
386 
387 
388 
390  if (m_output) cerr << "\rChecking for valid XML [68%]";
391  school_work_config = pop_config.get_child("school_work_profile.employable");
392  auto work_config = pop_config.get_child("work");
393 
394  total_size = work_config.get<int>("<xmlattr>.size");
395  minimum = school_work_config.get<int>("employee.<xmlattr>.min");
396  max = school_work_config.get<int>("employee.<xmlattr>.max");
397  fraction = school_work_config.get<double>("<xmlattr>.fraction") / 100.0;
398 
399  if (minimum > max || minimum < 0 || max < 0) {
400  throw invalid_argument("In PopulationGenerator: Numerical error in min/max pair.");
401  }
402 
403  if (total_size <= 0 || cluster_size <= 0 || fraction < 0.0 || fraction > 1.0) {
404  throw invalid_argument("In PopulationGenerator: Numerical error.");
405  }
406 
407  int min2 = school_work_config.get<int>("young_employee.<xmlattr>.min");
408  int max2 = school_work_config.get<int>("young_employee.<xmlattr>.max");
409  fraction = work_config.get<double>("far.<xmlattr>.fraction") / 100.0;
410 
411  if (min2 > max2 || min2 < 0 || max2 < 0) {
412  throw invalid_argument("In PopulationGenerator: Numerical error in min/max pair.");
413  }
414 
415  if (max2 >= minimum) {
416  throw invalid_argument("In PopulationGenerator: Overlapping min/max pairs.");
417  }
418 
419  if (fraction < 0.0 || fraction > 1.0) {
420  throw invalid_argument("In PopulationGenerator: Numerical error.");
421  }
422 
423 
424 
426  if (m_output) cerr << "\rChecking for valid XML [84%]";
427  total_size = pop_config.get<int>("community.<xmlattr>.size");
428 
429  if (total_size <= 0) {
430  throw invalid_argument("In PopulationGenerator: Numerical error.");
431  }
432  if (m_output) cerr << "\rChecking for valid XML [100%]\n";
433 }
434 
435 template<class U>
437  m_next_id = 1;
438  string file_name = m_props.get<string>("population.family.<xmlattr>.file");
439 
440  FamilyParser parser;
441 
442  vector<FamilyConfig> family_config {parser.parseFamilies(file_name)};
443 
444  uint current_generated = 0;
445 
447  AliasDistribution dist {vector<double>(family_config.size(), 1.0 / family_config.size())};
448  while (current_generated < m_total) {
449  if (m_output)
450  cerr << "\rGenerating households [" << min(uint(double(current_generated) / m_total * 100), 100U) << "%]";
451 
453  uint family_index = dist(m_rng);
454  FamilyConfig& new_config = family_config.at(family_index);
455 
457  SimpleHousehold new_household;
458  new_household.m_id = m_next_id;
459  for (uint& age: new_config) {
460  SimplePerson new_person {age, m_next_id};
461  new_person.m_household_id = m_next_id;
462  m_people.push_back(new_person);
463  new_household.m_indices.push_back(m_people.size() - 1);
464 
466  m_age_distribution[age] = m_age_distribution[age] + 1;
467  }
468 
470  m_household_size[new_config.size()] = m_household_size[new_config.size()] + 1;
471 
472  m_households.push_back(new_household);
473  m_next_id++;
474  current_generated += new_config.size();
475  }
476  if (m_output) cerr << "\rGenerating households [100%]...\n";
477 }
478 
479 template<class U>
481  ptree cities_config = m_props.get_child("population.cities");
482  m_next_id = 1;
483  uint size_check = 0;
484 
485  uint generated = 0;
486  uint to_generate = distance(cities_config.begin(), cities_config.end());
487  for (auto it = cities_config.begin(); it != cities_config.end(); it++) {
488  if (m_output) cerr << "\rGenerating cities [" << min(uint(double(generated) / to_generate * 100), 100U) << "%]";
489  if (it->first == "city") {
490  string name = it->second.get<string>("<xmlattr>.name");
491  int size = it->second.get<int>("<xmlattr>.pop");
492  double latitude = it->second.get<double>("<xmlattr>.lat");
493  double longitude = it->second.get<double>("<xmlattr>.lon");
494  size_check += size;
495 
496  SimpleCity new_city = SimpleCity(0, size, m_next_id, name, GeoCoordinate(latitude, longitude));
497  ++m_next_id;
498 
499  m_cities.push_back(new_city);
500  }
501  generated++;
502  }
503  if (m_output) cerr << "\rGenerating cities [100%]...\n";
504 
506  auto compare_city_size = [](const SimpleCity& a, const SimpleCity b) { return a.m_max_size > b.m_max_size; };
507  sort(m_cities.begin(), m_cities.end(), compare_city_size);
508 }
509 
510 template<class U>
512  double latitude_middle = 0.0;
513  double longitude_middle = 0.0;
514  for (const SimpleCity& city: m_cities) {
515  latitude_middle += city.m_coord.m_latitude;
516  longitude_middle += city.m_coord.m_longitude;
517  }
518  latitude_middle /= m_cities.size();
519  longitude_middle /= m_cities.size();
520 
521  GeoCoordinate result;
522  result.m_latitude = latitude_middle;
523  result.m_longitude = longitude_middle;
524  return result;
525 }
526 
527 template<class U>
529  double current_maximum = -1.0;
530 
532  for (const SimpleCity& city: m_cities) {
533  double distance = calc.getDistance(coord, city.m_coord);
534  if (distance > current_maximum) {
535  current_maximum = distance;
536  }
537  }
538 
539  return current_maximum;
540 }
541 
542 template<class U>
544  uint result = 0;
545 
546  for (const SimpleCity& city: m_cities) {
547  result += city.m_max_size;
548  }
549 
550  return result;
551 }
552 
553 template<class U>
555  uint result = 0;
556 
557  for (const SimpleCluster& village: m_villages) {
558  result += village.m_max_size;
559  }
560 
561  return result;
562 }
563 
564 template<class U>
566  // Do NOT reset the id counter (cities and villages will be treated as one)
567  ptree village_config = m_props.get_child("population.villages");
568  double village_radius_factor = village_config.get<double>("<xmlattr>.radius");
569  GeoCoordinate middle = getCityMiddle();
570  double radius = getCityRadius(middle);
571  uint city_population = getCityPopulation();
572  int unassigned_population = m_people.size() - city_population;
573  int unassigned_population_progress = unassigned_population;
574 
576  vector<double> fractions;
577  vector<MinMax> boundaries;
578  for (auto it = village_config.begin(); it != village_config.end(); it++) {
579  if (it->first == "village") {
580  uint min = it->second.get<uint>("<xmlattr>.min");
581  uint max = it->second.get<uint>("<xmlattr>.max");
582  double fraction = it->second.get<double>("<xmlattr>.fraction") / 100.0;
583  MinMax min_max {min, max};
584 
585  fractions.push_back(fraction);
586  boundaries.push_back(min_max);
587  }
588  }
589 
591  AliasDistribution village_type_dist {fractions};
593  while (unassigned_population > 0) {
594  if (m_output)
595  cerr << "\rGenerating villages [" << min(uint(
596  double(unassigned_population_progress - unassigned_population) / unassigned_population_progress *
597  100), 100U) << "%]";
598  uint village_type_index = village_type_dist(m_rng);
599  MinMax village_pop = boundaries.at(village_type_index);
600  uint range_interval_size = village_pop.max - village_pop.min + 1;
601 
602  AliasDistribution village_size_dist {vector<double>(range_interval_size, 1.0 / range_interval_size)};
603  uint village_size = village_size_dist(m_rng) + village_pop.min;
604 
605  SimpleCluster new_village;
606  new_village.m_max_size = village_size;
607  new_village.m_id = m_next_id;
608  m_next_id++;
609  new_village.m_coord = calc.generateRandomCoord(middle, radius * village_radius_factor, m_rng);
610 
612  auto same_coordinate_village = [&](const SimpleCluster& cl) { return cl.m_coord == new_village.m_coord; };
613  auto same_coordinate_city = [&](const SimpleCity& cl) { return cl.m_coord == new_village.m_coord; };
614 
615  auto it_villages = find_if(m_villages.begin(), m_villages.end(), same_coordinate_village);
616  auto it_cities = find_if(m_cities.begin(), m_cities.end(), same_coordinate_city);
617 
618  if (it_villages == m_villages.end() && it_cities == m_cities.end()) {
619  m_villages.push_back(new_village);
620  unassigned_population -= new_village.m_max_size;
621  }
622  }
623  if (m_output) cerr << "\rGenerating villages [100%]...\n";
624 }
625 
626 template<class U>
628  uint city_pop = getCityPopulation();
629  uint village_pop = getVillagePopulation();
630  uint total_pop =
631  village_pop + city_pop;
632 
635  vector<double> fractions;
636  for (SimpleCity& city: m_cities) {
637  fractions.push_back(double(city.m_max_size) / double(total_pop));
638  }
639 
640  for (SimpleCluster& village: m_villages) {
641  fractions.push_back(double(village.m_max_size) / double(total_pop));
642  }
643 
644  AliasDistribution village_city_dist {fractions};
645  int i = 0;
646  for (SimpleHousehold& household: m_households) {
647  if (m_output)
648  cerr << "\rPlacing households [" << min(uint(double(i) / m_households.size() * 100), 100U) << "%]";
649  uint index = village_city_dist(m_rng);
650  if (index < m_cities.size()) {
652  SimpleCity& city = m_cities.at(index);
653  city.m_current_size += household.m_indices.size();
654  for (uint& person_index: household.m_indices) {
655  m_people.at(person_index).m_coord = city.m_coord;
656  }
657  m_locations[make_pair(ClusterType::Household, household.m_id)] = city.m_coord;
658 
659  } else {
661  SimpleCluster& village = m_villages.at(index - m_cities.size());
662  village.m_current_size += household.m_indices.size();
663  for (uint& person_index: household.m_indices) {
664  m_people.at(person_index).m_coord = village.m_coord;
665  }
666  m_locations[make_pair(ClusterType::Household, household.m_id)] = village.m_coord;
667  }
668  i++;
669  }
670  if (m_output) cerr << "\rPlacing households [100%]...\n";
671 }
672 
673 template<class U>
676  m_next_id = 1;
677  ptree education_config = m_props.get_child("population.education");
678  ptree school_work_config = m_props.get_child("population.school_work_profile.mandatory");
679  uint school_size = education_config.get<uint>("mandatory.<xmlattr>.total_size");
680  uint min_age = school_work_config.get<uint>("<xmlattr>.min");
681  uint max_age = school_work_config.get<uint>("<xmlattr>.max");
682  uint cluster_size = education_config.get<uint>("mandatory.<xmlattr>.cluster_size");
683 
684  placeClusters(school_size, min_age, max_age, 1.0, m_mandatory_schools, "schools", ClusterType::School, false);
685 
686  // Split the schools in clusters
687  m_next_id = 1;
688  for (SimpleCluster& cluster: m_mandatory_schools) {
689  uint current_size = 0;
690 
691  m_mandatory_schools_clusters.push_back(vector<SimpleCluster>());
692  while (current_size < cluster.m_max_size) {
693  current_size += cluster_size;
694  SimpleCluster new_cluster;
695  new_cluster.m_max_size = cluster_size;
696  new_cluster.m_id = m_next_id;
697  m_next_id++;
698  new_cluster.m_coord = cluster.m_coord;
699  m_mandatory_schools_clusters.back().push_back(new_cluster);
700  m_locations[make_pair(ClusterType::School, new_cluster.m_id)] = new_cluster.m_coord;
701  }
702  }
703 }
704 
705 template<class U>
707  // Note: don't reset the next_id, the universities are also "schools", and the previously handled function was "makeSchools"
708  ptree school_work_config = m_props.get_child("population.school_work_profile.employable.young_employee");
709  ptree university_config = m_props.get_child("population.education.optional");
710  uint min_age = school_work_config.get<uint>("<xmlattr>.min");
711  uint max_age = school_work_config.get<uint>("<xmlattr>.max");
712  double fraction = 1.0 - school_work_config.get<double>("<xmlattr>.fraction") / 100.0;
713  uint size = university_config.get<uint>("<xmlattr>.total_size");
714  uint cluster_size = university_config.get<uint>("<xmlattr>.cluster_size");
715 
716  uint intellectual_pop = 0;
717  for (uint i = min_age; i <= max_age; i++) {
718  intellectual_pop += m_age_distribution[i];
719  }
720 
721  intellectual_pop = ceil(intellectual_pop * fraction);
722 
723  uint needed_universities = ceil(double(intellectual_pop) / size);
724  uint placed_universities = 0;
725  uint clusters_per_univ = size / cluster_size;
726  uint left_over_cluster_size = size % cluster_size;
727 
728  m_optional_schools.clear();
729 
730  while (needed_universities > placed_universities) {
731  if (m_output)
732  cerr << "\rPlacing Universities ["
733  << min(uint(double(needed_universities) / placed_universities * 100), 100U) << "%]";
734 
737  vector<SimpleCluster> univ;
738  for (uint i = 0; i < clusters_per_univ; i++) {
739  SimpleCluster univ_cluster;
740  univ_cluster.m_id = m_next_id;
741  univ_cluster.m_max_size = cluster_size;
742  univ_cluster.m_coord = m_cities.at(placed_universities % m_cities.size()).m_coord;
743  m_next_id++;
744  univ.push_back(univ_cluster);
745 
746  m_locations[make_pair(ClusterType::School, univ_cluster.m_id)] = univ_cluster.m_coord;
747  }
748 
749  if (left_over_cluster_size > 0) {
751  SimpleCluster univ_cluster;
752  univ_cluster.m_id = m_next_id;
753  univ_cluster.m_max_size = left_over_cluster_size;
754  univ_cluster.m_coord = m_cities.at(placed_universities % m_cities.size()).m_coord;
755  m_next_id++;
756  univ.push_back(univ_cluster);
757 
758  m_locations[make_pair(ClusterType::School, univ_cluster.m_id)] = univ_cluster.m_coord;
759  }
760 
761  m_optional_schools.push_back(univ);
762  placed_universities++;
763  }
764 }
765 
766 template<class U>
769  vector<SimpleCluster> result;
770 
771  for (SimpleCity& city: m_cities) {
772  for (SimpleCluster& workplace: m_workplaces) {
773  if (city.m_coord == workplace.m_coord) {
774  result.push_back(workplace);
775  }
776  }
777  }
778 
779  for (SimpleCluster& village: m_villages) {
780  for (SimpleCluster& workplace: m_workplaces) {
781  if (village.m_coord == workplace.m_coord) {
782  result.push_back(workplace);
783  }
784  }
785  }
786 
787  m_workplaces = result;
788 }
789 
790 template<class U>
792  m_next_id = 1;
793  ptree school_work_config = m_props.get_child("population.school_work_profile.employable");
794  ptree work_config = m_props.get_child("population.work");
795 
796  uint size = work_config.get<uint>("<xmlattr>.size");
797  uint min_age = school_work_config.get<uint>("employee.<xmlattr>.min");
798  uint max_age = school_work_config.get<uint>("employee.<xmlattr>.max");
799  uint young_min_age = school_work_config.get<uint>("young_employee.<xmlattr>.min");
800  uint young_max_age = school_work_config.get<uint>("young_employee.<xmlattr>.max");
801  double young_fraction = school_work_config.get<double>("young_employee.<xmlattr>.fraction") / 100.0;
802  double fraction = school_work_config.get<double>("<xmlattr>.fraction") / 100.0;
803 
804  uint possible_students = 0;
805  uint total_old = 0;
806  for (SimplePerson& person: m_people) {
807  if (person.m_age >= young_min_age && person.m_age <= young_max_age) {
808  possible_students++;
809  }
810 
811  if (person.m_age >= std::max(min_age, young_max_age + 1) && person.m_age <= max_age) {
812  total_old++;
813  }
814  }
815 
816  uint total_of_age = possible_students + total_old;
817  // Subtract actual students
818  uint total_working = total_of_age - ceil(possible_students * (1.0 - young_fraction));
819  total_working = ceil(total_working * fraction);
820 
821  // Calculate the actual fraction of people between young_min_age and max_age who are working
822  double actual_fraction = double(total_working) / total_of_age;
823 
824  placeClusters(size, young_min_age, max_age, actual_fraction, m_workplaces, "workplaces", ClusterType::Work);
825 
827  sortWorkplaces();
828 }
829 
830 template<class U>
833  m_next_id = 1;
834  ptree community_config = m_props.get_child("population.community");
835  uint size = community_config.get<uint>("<xmlattr>.size");
836 
837  placeClusters(size, 0, 0, 1.0, m_primary_communities, "primary communities", ClusterType::PrimaryCommunity);
838  m_next_id = 1;
839  placeClusters(size, 0, 0, 1.0, m_secondary_communities, "secondary communities", ClusterType::SecondaryCommunity);
840 }
841 
842 template<class U>
844  m_next_id = 1;
845  ptree education_config = m_props.get_child("population.education.mandatory");
846  ptree school_work_config = m_props.get_child("population.school_work_profile.mandatory");
847  uint min_age = school_work_config.get<uint>("<xmlattr>.min");
848  uint max_age = school_work_config.get<uint>("<xmlattr>.max");
849  double start_radius = m_props.get<double>("population.commutingdata.<xmlattr>.start_radius");
850 
851  double factor = m_props.get<double>("population.commutingdata.<xmlattr>.factor");
852 
853  double current_radius = start_radius;
854 
855  uint total = 0;
856  uint total_placed = 0;
857 
858  for (uint i = min_age; i <= max_age; i++) {
859  total += m_age_distribution[i];
860  }
861 
862  auto distance_map = makeDistanceMap(start_radius, factor, m_mandatory_schools);
863 
864  for (SimplePerson& person: m_people) {
865  current_radius = start_radius;
866  if (person.m_age >= min_age && person.m_age <= max_age) {
867  if (m_output)
868  cerr << "\rAssigning children to schools [" << min(uint(double(total_placed) / total * 100), 100U)
869  << "%]";
870  total_placed++;
871 
872  vector<uint> closest_clusters_indices;
873 
874  while (closest_clusters_indices.size() == 0 && m_mandatory_schools.size() != 0) {
875  closest_clusters_indices = getClustersWithinRange(current_radius, distance_map, person.m_coord);
876  current_radius *= factor;
877  }
878 
879  AliasDistribution cluster_dist {
880  vector<double>(closest_clusters_indices.size(), 1.0 / double(closest_clusters_indices.size()))};
881  uint index = closest_clusters_indices.at(cluster_dist(m_rng));
882 
883  AliasDistribution inner_cluster_dist {vector<double>(m_mandatory_schools_clusters.at(index).size(),
884  1.0 / m_mandatory_schools_clusters.at(index).size())};
885  uint index2 = inner_cluster_dist(m_rng);
886 
887  m_mandatory_schools_clusters.at(index).at(index2).m_current_size++;
888  person.m_school_id = m_mandatory_schools_clusters.at(index).at(index2).m_id;
889  }
890  }
891  if (m_output) cerr << "\rAssigning children to schools [100%]...\n";
892 }
893 
894 template<class U>
896  ptree school_work_config = m_props.get_child("population.school_work_profile.employable.young_employee");
897  ptree university_config = m_props.get_child("population.education.optional");
898  uint min_age = school_work_config.get<uint>("<xmlattr>.min");
899  uint max_age = school_work_config.get<uint>("<xmlattr>.max");
900  double student_fraction = 1.0 - school_work_config.get<double>("<xmlattr>.fraction") / 100.0;
901  double commute_fraction = university_config.get<double>("far.<xmlattr>.fraction") / 100.0;
902  double radius = m_props.get<double>("population.commutingdata.<xmlattr>.start_radius");
903 
904  AliasDistribution commute_dist {{commute_fraction, 1.0 - commute_fraction}};
905  AliasDistribution student_dist {{student_fraction, 1.0 - student_fraction}};
906 
907  uint total = 0;
908  uint total_placed = 0;
909 
910  for (uint i = min_age; i <= max_age; i++) {
911  total += m_age_distribution[i];
912  }
913 
914  auto distance_map = makeDistanceMap(radius, 2.0, m_optional_schools);
915 
916  for (SimplePerson& person: m_people) {
917  if (person.m_age >= min_age && person.m_age <= max_age && student_dist(m_rng) == 0) {
918  if (m_output)
919  cerr << "\rAssigning students to universities [" << min(uint(double(total_placed) / total * 100), 100U)
920  << "%]";
921  total_placed++;
922 
923  if (commute_dist(m_rng) == 0) {
925  assignCommutingStudent(person, distance_map);
926  } else {
928  assignCloseStudent(person, radius, distance_map);
929  }
930  }
931  }
932  if (m_output) cerr << "\rAssigning students to universities [100%]...\n";
933 }
934 
935 template<class U>
936 void PopulationGenerator<U>::removeFromUniMap(vector<pair<GeoCoordinate, map<double, vector<uint>>>>& distance_map,
937  uint index) const {
938  bool full_capacity = true;
939 
940  for (uint curr_univ = index; curr_univ < m_optional_schools.size(); ++curr_univ) {
941  for (const SimpleCluster& uni: m_optional_schools.at(curr_univ)) {
942  if (uni.m_current_size < uni.m_max_size) {
943  full_capacity = false;
944  }
945  }
946  }
947 
949  if (full_capacity) {
950  for (auto& coord_map_pair: distance_map) {
951  for (auto it = coord_map_pair.second.begin(); it != coord_map_pair.second.end(); ++it) {
952  auto& index_vector = it->second;
953 
954  auto cluster_iterator = std::find(index_vector.begin(), index_vector.end(), index);
955 
956  if (cluster_iterator != index_vector.end()) {
957  index_vector.erase(cluster_iterator, cluster_iterator + 1);
958 
959  removeFromUniMap(distance_map, index);
960 
961  return;
962  }
963 
964  }
965  }
966  }
967 }
968 
969 template<class U>
970 bool PopulationGenerator<U>::removeFromMap(vector<pair<GeoCoordinate, map<double, vector<uint>>>>& distance_map,
971  uint index) const {
972 
973  // Remove this cluster from the distance map
974  for (auto& coord_map_pair: distance_map) {
975  for (auto it = coord_map_pair.second.begin(); it != coord_map_pair.second.end(); ++it) {
976  auto& index_vector = it->second;
977 
978  auto cluster_iterator = std::find(index_vector.begin(), index_vector.end(), index);
979 
980  if (cluster_iterator != index_vector.end()) {
981  index_vector.erase(cluster_iterator, cluster_iterator + 1);
982 
983  removeFromMap(distance_map, index);
984  return true;
985  }
986  }
987  }
988 
989  return false;
990 }
991 
992 template<class U>
994  vector<pair<GeoCoordinate, map<double, vector<uint>>>>& distance_map) {
995  uint current_city = 0;
996  bool added = false;
997 
998  while (current_city < m_cities.size() && !added) {
999  uint current_univ = current_city;
1000  while (current_univ < m_optional_schools.size() && !added) {
1001  for (uint i = 0; i < m_optional_schools.at(current_univ).size(); ++i) {
1002  SimpleCluster& univ_cluster = m_optional_schools.at(current_univ).at(i);
1003  if (univ_cluster.m_current_size < univ_cluster.m_max_size) {
1004  univ_cluster.m_current_size++;
1005  person.m_school_id = univ_cluster.m_id;
1006  added = true;
1007 
1008  removeFromUniMap(distance_map, current_city);
1009  break;
1010  }
1011  }
1012  current_univ += m_cities.size();
1013  }
1014  current_city++;
1015  }
1016 }
1017 
1018 template<class U>
1020  vector<pair<GeoCoordinate, map<double, vector<uint>>>>& distance_map) {
1021  double factor = m_props.get<double>("population.commutingdata.<xmlattr>.factor");
1022  double current_radius = start_radius;
1023  bool added = false;
1024  vector<uint> closest_clusters_indices;
1025 
1026  while (!added) {
1027 
1028  // Note how getting the distance to the closest univ is the same as getting the distance to the closest city
1029  // This is because their vectors are in the same order!
1030 
1031  closest_clusters_indices = getClustersWithinRange(current_radius, distance_map, person.m_coord);
1032 
1033  if (closest_clusters_indices.size() != 0) {
1034  AliasDistribution dist {
1035  vector<double>(closest_clusters_indices.size(), 1.0 / closest_clusters_indices.size())};
1036 
1037  uint index = dist(m_rng);
1038  while (index < m_optional_schools.size() && !added) {
1039 
1040  for (uint i = 0; i < m_optional_schools.at(index).size(); i++) {
1041  SimpleCluster& univ_cluster = m_optional_schools.at(index).at(i);
1042  if (univ_cluster.m_current_size < univ_cluster.m_max_size) {
1043  univ_cluster.m_current_size++;
1044  person.m_school_id = univ_cluster.m_id;
1045  added = true;
1046 
1047  removeFromUniMap(distance_map, index);
1048  break;
1049  }
1050  }
1051  index += m_cities.size();
1052  }
1053 
1054  if (added) {
1055  break;
1056  }
1057  }
1058 
1059  current_radius *= factor;
1060  }
1061 }
1062 
1063 template<class U>
1065  ptree school_work_config = m_props.get_child("population.school_work_profile.employable");
1066  ptree work_config = m_props.get_child("population.work");
1067  uint min_age = school_work_config.get<uint>("young_employee.<xmlattr>.min");
1068  uint max_age = school_work_config.get<uint>("employee.<xmlattr>.max");
1069  double unemployment_rate = 1.0 - school_work_config.get<double>("<xmlattr>.fraction") / 100.0;
1070  double commute_fraction = work_config.get<double>("far.<xmlattr>.fraction") / 100.0;
1071  double radius = m_props.get<double>("population.commutingdata.<xmlattr>.start_radius");
1072 
1073  AliasDistribution unemployment_dist {{unemployment_rate, 1.0 - unemployment_rate}};
1074  AliasDistribution commute_dist {{commute_fraction, 1.0 - commute_fraction}};
1075 
1076  uint total = 0;
1077  uint total_placed = 0;
1078  uint amount_full = 0;
1079 
1080  for (uint i = min_age; i <= max_age; i++) {
1081  total += m_age_distribution[i];
1082  }
1083 
1084  auto distance_map = makeDistanceMap(radius, 2.0, m_workplaces);
1085 
1086  for (SimplePerson& person: m_people) {
1087  if (amount_full == m_workplaces.size()) {
1088  // No more workplaces left, just stop
1089  break;
1090  }
1091 
1092  if (person.m_age >= min_age && person.m_age <= max_age) {
1093  if (m_output)
1094  cerr << "\rAssigning people to workplaces [" << min(uint(double(total_placed) / total * 100), 100U)
1095  << "%]";
1096  total_placed++;
1097  if (unemployment_dist(m_rng) == 1 && person.m_school_id == 0) {
1098  bool filled_cluster = true;
1099  if (commute_dist(m_rng) == 0) {
1101  filled_cluster = assignCommutingEmployee(person, distance_map);
1102  } else {
1104  filled_cluster = assignCloseEmployee(person, radius, distance_map);
1105  }
1106 
1107  if (filled_cluster) {
1108  ++amount_full;
1109  }
1110  }
1111  }
1112  }
1113  if (m_output) cerr << "\rAssigning people to workplaces [100%]...\n";
1114 }
1115 
1116 template<class U>
1118  vector<pair<GeoCoordinate, map<double, vector<uint>>>>& distance_map) {
1122  for (uint i = 0; i < m_workplaces.size(); ++i) {
1123  SimpleCluster& workplace = m_workplaces.at(i);
1124 
1125  if (workplace.m_max_size > workplace.m_current_size) {
1126  workplace.m_current_size++;
1127  person.m_work_id = workplace.m_id;
1128 
1129  if (workplace.m_current_size >= workplace.m_max_size) {
1130  return removeFromMap(distance_map, i);
1131  }
1132 
1133  break;
1134  }
1135  }
1136 
1137  return false;
1138 }
1139 
1140 template<class U>
1142  vector<pair<GeoCoordinate, map<double, vector<uint>>>>& distance_map) {
1143  double factor = m_props.get<double>("population.commutingdata.<xmlattr>.factor");
1144  double current_radius = start_radius;
1145 
1146  vector<uint> closest_clusters_indices;
1147 
1148  while (closest_clusters_indices.size() == 0) {
1149  closest_clusters_indices = getClustersWithinRange(current_radius, distance_map, person.m_coord);
1150  current_radius *= factor;
1151  }
1152 
1153  AliasDistribution dist {
1154  vector<double>(closest_clusters_indices.size(), 1.0 / double(closest_clusters_indices.size()))};
1155  uint rnd = dist(m_rng);
1156  auto index = closest_clusters_indices.at(rnd);
1157  SimpleCluster& workplace = m_workplaces.at(index);
1158 
1159  person.m_work_id = workplace.m_id;
1160 
1161  workplace.m_current_size++;
1162 
1163  if (workplace.m_current_size >= workplace.m_max_size) {
1164  return removeFromMap(distance_map, index);
1165  }
1166 
1167 
1168  return false;
1169 }
1170 
1171 template<class U>
1172 void PopulationGenerator<U>::assignToCommunities(vector<pair<GeoCoordinate, map<double, vector<uint>>>>& distance_map,
1173  vector<SimpleCluster>& clusters,
1174  uint SimplePerson::* member,
1175  const string& name) {
1176 
1177  double start_radius = m_props.get<double>("population.commutingdata.<xmlattr>.start_radius");
1178  double factor = m_props.get<double>("population.commutingdata.<xmlattr>.factor");
1179 
1180  uint total = m_households.size();
1181  uint total_placed = 0;
1182 
1183  for (SimpleHousehold& household: m_households) {
1184  if (m_output)
1185  cerr << "\rAssigning people to " << name << " [" << min(uint(double(total_placed) / total * 100), 100U)
1186  << "%]";
1187  total_placed++;
1188 
1189  double current_radius = start_radius;
1190  vector<uint> closest_clusters_indices;
1191 
1192  // TODO make sure merge didn't screw this up
1193  while (closest_clusters_indices.size() == 0) {
1194  closest_clusters_indices = getClustersWithinRange(current_radius, distance_map,
1195  m_people.at(household.m_indices.back()).m_coord);
1196  current_radius *= factor;
1197  }
1198 
1199  AliasDistribution dist {
1200  vector<double>(closest_clusters_indices.size(), 1.0 / double(closest_clusters_indices.size()))};
1201  uint index = closest_clusters_indices.at(dist(m_rng));
1202  SimpleCluster& community = clusters.at(index);
1203  for (uint& person_index: household.m_indices) {
1204  SimplePerson& person = m_people.at(person_index);
1205  person.*member = community.m_id;
1206  community.m_current_size++;
1207  }
1208 
1210  if (community.m_current_size >= community.m_max_size) {
1211  removeFromMap(distance_map, index);
1212  }
1213 
1214  }
1215  if (m_output) cerr << "\rAssigning people to " << name << " [100%]...\n";
1216 }
void writeClusters(const string &target_clusters) const
Writes the clusters to the file (type, ID and coordinates), see PopulationGenerator::generate.
void assignCommutingStudent(SimplePerson &person, vector< pair< GeoCoordinate, map< double, vector< uint >>>> &distance_map)
Put one student in a university according to the rules of commuting students.
Interface for install directory queries.
string toString(ClusterType c)
Converts a ClusterType value to corresponding name.
Definition: ClusterType.cpp:54
bool removeFromMap(vector< pair< GeoCoordinate, map< double, vector< uint >>>> &distance_map, uint index) const
Remove an element from the map (of regular clusters, not like universities, because they represent a ...
bool assignCloseEmployee(SimplePerson &person, double start_radius, vector< pair< GeoCoordinate, map< double, vector< uint >>>> &distance_map)
Assign one person to a workplace according to the rule of workers that work close to their home...
unsigned int uint
Definition: FamilyParser.h:12
Usage is very simple, construct with a vector of probabilities, then use as a distribution from the s...
void assignToUniversities()
Put students in universities.
void writePop(const string &target_pop) const
Writes the population to the file, see PopulationGenerator::generate.
double getVillagePopulation() const
Get the population of all the villages combined (if they were on full capacity)
Time Dependent Person DataType.
Definition: NoBehaviour.h:17
GeoCoordinate generateRandomCoord(const GeoCoordinate &coord, double radius, T rng) const
Result is in kilometers Uses the haversine formula See: http://www.movable-type.co.uk/scripts/latlong.html.
static boost::filesystem::path getDataDir()
Utility method: get path to the directory for data files.
vector< FamilyConfig > parseFamilies(string filename) const
void assignToWork()
Assign people to a workplace.
void removeFromUniMap(vector< pair< GeoCoordinate, map< double, vector< uint >>>> &distance_map, uint index) const
Remove an element from the university map (the university map is special compared to other cluster ma...
double getCityRadius(const GeoCoordinate &coord) const
Returns the distance between the given coordinate and the furthest city (in km)
util::GeoCoordinate m_coord
Definition: popgen/utils.h:50
void generate(const string &prefix)
Generates a population, writes the result to the files found in the data directory Output files are r...
void writeCities(const string &target_cities)
Writes the cities to the file, see PopulationGenerator::generate, recently, the villages have been ad...
void makeWork()
Make workplaces Note: due to the fact that the cluster size of workplaces must be respected...
util::GeoCoordinate m_coord
Definition: popgen/utils.h:67
void makeUniversities()
Make the universities, place them in a city.
void checkForValidXML() const
Checks the xml on correctness, this includes only semantic errors, no syntax errors.
std::vector< uint > m_indices
Definition: popgen/utils.h:43
bool assignCommutingEmployee(SimplePerson &person, vector< pair< GeoCoordinate, map< double, vector< uint >>>> &distance_map)
Assign one person to a workplace according to the rule of commuting workers.
void writeHouseholds(const string &target_households) const
Writes the households to the file, see PopulationGenerator::generate.
void assignToCommunities(vector< pair< GeoCoordinate, map< double, vector< uint >>>> &distance_map, vector< SimpleCluster > &clusters, uint SimplePerson::*member, const string &name="")
Assign entire households.
void makeSchools()
Make the schools, place them in a village/city.
void makeHouseholds()
Generates all households (not yet their positions)
void makeCommunities()
Make the communities.
PopulationGenerator(const string &filename, const int &seed, bool output=true)
Constructor: Check if the xml is valid and set up the basic things like a random generator.
vector< uint > FamilyConfig
Definition: FamilyParser.h:13
void assignCloseStudent(SimplePerson &person, double start_radius, vector< pair< GeoCoordinate, map< double, vector< uint >>>> &distance_map)
Put one student in a university according to the rules of students that study close to their home...
void sortWorkplaces()
Make sure workplaces are sorted: workplaces in bigger cities are in front of workplaces in smaller ci...
util::GeoCoordinate m_coord
Definition: popgen/utils.h:37
void makeVillages()
Generate all villages (without inhabitants)
void makeCities()
Generate all cities (without inhabitants)
double getDistance(const GeoCoordinate &coord1, const GeoCoordinate &coord2) const
static const GeoCoordCalculator & getInstance()
Singleton pattern.
void assignToSchools()
Put children in mandatory schools.
TimeStamp class.
GeoCoordinate getCityMiddle() const
Gets the middle of all cities.
double getCityPopulation() const
Get the population of all the cities combined (if they were on full capacity)
void placeHouseholds()
Assign the households to a city/village.