I was curious a few weeks ago about how surnames come and go (mostly on the go). I wanted to know how quickly surnames die out, because (in my simplified simulation world) once the last person with a surname dies, it's gone forever. I started writing a small simulation to see this happen in "real time". I started with a simple "person" representation that behaved exactly like all the other "persons" in the simulation. However, as I got into the project, I started coming up with certain traits, since not everyone lives to the same age, not everyone has the same number of children, not everyone pairs up to even have children, etc.
A few assumptions:
- I didn't base this simulation on any actual population statistics, though I think it does a decent job at producing some interesting results.
- Surnames do not regenerate. The starting population each gets a unique "name" (an integer one greater than the last), and that forms the maximum amount of names in the simulation.
- This is paternally biased. To keep things simple, females of the population take the male surname when paired. Surnames are never combined (see point 2), or modified in any way. This simulation would work the same (mostly?) if maternally biased, though I suppose at that point males are only needed for population growth (topic not relevant for this simulation, just some side thoughts).
- Surnames are inherited from the parent. I didn't originally consider this an assumption, but according to Toby Speight, this is not always true in every culture, with Iceland being the example.
- People die when they reach their maximum lifespan, but can also die early due to "accidents", being more prone the old they are.
- When a person with a partner dies, the surviving partner does not gain a new partner.
- Exponential distributions are used here to help move away from uniform distributions, since I found that uniform populations die really easily. An example would be that a
Personis more likely to have a life span of 80 than 20, but that drops off significantly past 80.
Therefore, I came up with the following
Person object to act as the individual units of my simulation:
# Person.py import random import numpy as np male = 1 female = 2 # these reduce performance contributions from lookups by half, and because of how often they are called, any bit helps # # $ python -m timeit -s "import numpy as np" "np.random.randint(0, 1000)" # # 10000000 loops, best of 3: 0.168 usec per loop # # # $ python -m timeit -s "import numpy as np; npRandInt = np.random.randint" "npRandInt(0, 1000)" # # 10000000 loops, best of 3: 0.0814 usec per loop npRandInt = np.random.randint npRandUni = np.random.uniform npRandEx = np.random.exponential class Person(object): currentName = 0 deaths = 0 maxAgeDistribution =  maxChildDistribution = [0, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5] def __init__(self, parent=None): if parent is None: Person.currentName += 1 self.surname = Person.currentName else: self.surname = parent.surname self.gender = random.choice([male, female]) self.currentChildCount = 0 self.currentAge = 0 self.maxChildren = random.choice(Person.maxChildDistribution) # cap how many children one person can have self.maxAge = random.choice(Person.maxAgeDistribution) self.childDifficulty = npRandInt(0, 100 + 1) / 100.0 # how difficult is it to have a child self.childYearsStart = npRandInt(10, 15 + 1) self.childYearsEnd = npRandInt(30, 50 + 1) self.partner = None # some other person object self.withChild = False self.hasHadPartner = False # give a way to check if they have had a partner before, so if one dies, they don't get another # determine age person is eligible to gain a partner. Having this too high (on average) leads to excessive population death if self.childYearsStart < self.maxAge: self.partnerEligibleAge = self.childYearsStart + int(self.childYearsStart * npRandEx(0.2)) if self.partnerEligibleAge > self.maxAge: self.partnerEligibleAge = self.maxAge else: self.partnerEligibleAge = self.maxAge def __repr__(self): return 'Person(Gender: %6s, Surname: %5d, Age: %2d, MaxAge: %2d, Children: %1d, MaxChildren: %1d, Partner: %5s, PartnerEligibleAge: %2d)' % \ ("Male" if self.isMale() else "Female", self.surname, self.currentAge, self.maxAge, self.currentChildCount, self.maxChildren, "True" if self.partner is not None else "False", self.partnerEligibleAge) def gainPartner(self, partner): # gain a partner if both persons are eligible, female taking the surname of the other partner if self.isPartnerEligible(): if self.gender == female: self.surname = partner.surname self.partner = partner self.hasHadPartner = True def age(self): self.currentAge += 1 if self.currentAge >= self.maxAge: # if they outlive their lifespan, they die return self.die() percentAge = 1.0 - ((self.maxAge - self.currentAge) / float(self.maxAge)) # get their current age as a percent of their max age # get a random [0, 1] with an exponential distribution. # Magic numbers were determined experimentally, seeing what allowed population growth, and then by checking the distribution using pyplot randDeathPercent = 1.0 - npRandEx(0.05) * 2.0 # random accidents happen, so even if they haven't reached their max age, # there is some risk of death, especially as they get closer to their max age return self.die() if percentAge > randDeathPercent else False # perform self.die, which returns True, if they have died def gainChild(self): # gain a child if eligible and passes a "difficulty" check if self.isChildEligible(): chance = 1.0 - (self.childDifficulty * self.partner.childDifficulty) if chance > npRandUni(0, 1): self.withChild = True def haveChild(self, modifier=2.0): # produce child if they currently have one and pass a modified "difficulty" check if self.withChild is True: self.withChild = False # pass or fail, they will no longer have a child chance = 1.0 - self.childDifficulty if chance * modifier > npRandUni(0, 1): # if this modifier is set too low, not enough children are created and the population dies self.currentChildCount += 1 return Person(parent=self) # return a new person, giving it the last name of the parent return None # if we get this far, a child fails to be created def die(self): # if they have a partner, remove themselves from that person if self.partner is not None: self.partner.partner = None Person.deaths += 1 # add to the death stats return True def hasPartner(self): return self.partner is not None def isChildEligible(self): # check if person can have a child # Criteria: # only females can have children # they must have a partner # they haven't had more than their max children # they are within child bearing years, and so is their partner # they don't currently have a child return self.gender == female and \ self.partner is not None and \ self.currentChildCount < self.maxChildren and \ self.childYearsStart <= self.currentAge <= self.childYearsEnd and \ self.partner.childYearsStart <= self.partner.currentAge <= self.partner.childYearsEnd and \ self.withChild is False def isPartnerEligible(self): # check if person can gain a partner # Criteria: # they don't currently have a partner # they are old enough to have a partner # they haven't had a partner in the past return self.partner is None and \ self.currentAge >= self.partnerEligibleAge and \ self.hasHadPartner is False # helper functions to check for gender. # Performance profiling indicates these calls are expensive over time, so it is better to check this directly def isMale(self): return self.gender == male def isFemale(self): return self.gender == female
Now I needed to actually run the simulation. After a lot of performance checking, both profiling and seeing the results of the simulation, I was able to come up with the following main application, with
main handling the bulk of the administrative work, and
sim doing the actual simulation:
# main.py import random import numpy as np import cProfile import operator import datetime import math import matplotlib.pyplot as plt import cPickle as pickle import person from person import Person from stats import Stats import sys import uuid import time npRandInt = np.random.randint npRandUni = np.random.uniform npRandEx = np.random.exponential fileId = str(uuid.uuid1()) SERIALIZE = False def getTimeDelta(start, stop): delta = stop - start return (delta.microseconds + (delta.seconds + delta.days * 24 * 3600) * 10 ** 6) / (10 ** 6 * 1.0) * 1000.0 def setupAgeDistribution(distribution): for i in range(10000): distribution.append(int(getRandomAgeTo80())) for i in range(2000): distribution.append(int(getRandomAgeFrom80())) # get a random age in range [0, 80] # Magic numbers were determined experimentally, # seeing what allowed population growth, # and then by checking the distribution using pyplot def getRandomAgeTo80(): ran = 80.0 - npRandEx(50.0) * .4 while ran < 0: ran = 80.0 - npRandEx(50.0) * .4 return ran # get a random age in range [80, 100] # Magic numbers were determined experimentally, # seeing what allowed population growth, # and then by checking the distribution using pyplot def getRandomAgeFrom80(): ran = 80 + npRandEx(50.0) * .1 while ran > 100.0: ran = 80 + npRandEx(50.0) * .1 return ran def countItems(aList): # no access to collections.Counter, so use this instead array = np.bincount(np.array(aList)) nonZero = np.nonzero(array) return dict(zip(nonZero, array[nonZero])) def gainChildren(people): # have each person try to gain a child, letting the eligibility be determined by the object for p in people: p.gainChild() def gainPartners(people): # get a list of eligible partners, split by gender mPartners = [x for x in people if x.gender == person.male and x.isPartnerEligible()] fPartners = [x for x in people if x.gender == person.female and x.isPartnerEligible()] numPartners = min([len(mPartners), len(fPartners)]) # get the smaller count of partners if numPartners > 0: numChoices = npRandInt(0, numPartners + 1) # select a number of them to gain a partner if numChoices > 0: # shuffle both lists, splice to make lists the same length, zip them together, and pair up the zipped pairs. # Much faster than using random.choice. random.shuffle(mPartners) random.shuffle(fPartners) pairs = zip(mPartners[:numChoices], fPartners[:numChoices]) for m, f in pairs: m.gainPartner(f) f.gainPartner(m) def haveChildren(people): newCount = 0 for p in people.copy(): if p.withChild is True: newChild = p.haveChild(modifier=1.0) if newChild is not None: people.add(newChild) newCount += 1 return newCount def ageSim(people): for p in people.copy(): if p.age(): # if they ded people.remove(p) def sim(people): ageSim(people) # age the whole sim, removing those that die gainPartners(people) # pair up eligible people objects newCount = haveChildren(people) # have children before gaining them to "develop" them for a year gainChildren(people) # gain new children to attempt being created next year return newCount # return how many total new people have been added def main(people=None, stats=None): # arguments allow resuming from serialized (saved) data setupAgeDistribution(Person.maxAgeDistribution) random.shuffle(Person.maxAgeDistribution) random.shuffle(Person.maxChildDistribution) if people is None and stats is None: # set up the stats object, setting the initial population count, and how often to print the stats stats = Stats(oldDate=datetime.datetime.now(), startingPopulation=10000, statYear=100.0) people = set([Person() for i in range(stats.startingPopulation)]) # create a new set of random people stats.oldDate = datetime.datetime.now() try: for year in range(stats.year, stats.year + 500 + 1): stats.year = year if len(people) == 0: print "Population has died! Final stats below:\n" printStats(people, stats) break if year % stats.statYear == 0: printStats(people, stats) stats.newPersonCount += sim(people) except KeyboardInterrupt: # capture any ctrl-c, and print stats on exit printStats(people, stats) if len(people) > 0: sortedList = sorted(countItems([p.surname for p in people]).items(), key=operator.itemgetter(1)) # sort list by surname popularity if len(sortedList) < 1000: # only print it if the sorted list isn't too long. Console has an issue if the line is too long. sortedList.reverse() # display the list in descending order print sortedList sys.exit(0) # plot(people, stats) def printStats(people, stats): statYear = stats.statYear if stats.year % statYear != 0: statYear = stats.year % stats.statYear delta = getTimeDelta(stats.oldDate, datetime.datetime.now()) timeTaken = math.floor(delta / 1000.0) if timeTaken < 1: time.sleep(1) # don't spam the console output stats.timeTaken.append(timeTaken) population = len(people) stats.population.append(population) meanAge = math.floor(np.mean([person.currentAge for person in people])) if len(people) > 0 else 0 stats.meanAge.append(meanAge) stats.newPeople.append(stats.newPersonCount) peopleDelta = stats.newPersonCount - stats.oldPeople stats.oldPeople = stats.newPersonCount stats.deltaNewPeople.append(peopleDelta) surnameCount = len(countItems([p.surname for p in people])) if len(people) > 0 else 0 stats.surnames.append(surnameCount) stats.deaths.append(Person.deaths) deathsDelta = Person.deaths - stats.oldDeaths stats.oldDeaths = Person.deaths stats.deltaDeaths.append(deathsDelta) netPop = (peopleDelta / statYear) - (deathsDelta / statYear) stats.netPopDelta.append(netPop) pairCount = len([True for person in people if person.partner is not None]) if len(people) > 0 else 0 print datetime.datetime.now() print "Time taken: %d seconds (%d milliseconds / year)" % (timeTaken, math.floor(delta / statYear)) print "Current population size at year %d: %d" % (stats.year, population) print "Current mean age: %d" % meanAge print "Current new people: %d" % stats.newPersonCount print "Current deaths: %d" % Person.deaths print "Current delta new people: %d (%d people / year)" % (peopleDelta, (peopleDelta / statYear)) print "Current delta deaths: %d (%d deaths / year)" % (deathsDelta, (deathsDelta / statYear)) print "Current surname count: %d" % surnameCount print "Current partner percentage: %d" % (pairCount / float(population) * 100) if len(people) > 0 else 0 print "Net population per year: %d" % netPop print "" stats.oldDate = datetime.datetime.now() if SERIALIZE: # controls whether serialization is enabled, useful to turn off for when debugging or testing serialize(people, stats, fileId) def serialize(people, stats, id=""): with open("people%s.dat" % id, "wb") as f: pickle.dump(people, f, pickle.HIGHEST_PROTOCOL) with open("stats%s.dat" % id, "wb") as f: pickle.dump(stats, f, pickle.HIGHEST_PROTOCOL) def deserialize(people=True, stats=True, id=""): peopleData = None statsData = None if people: with open("people%s.dat" % id, "rb") as f: peopleData = pickle.load(f) if stats: with open("stats%s.dat" % id, "rb") as f: statsData = pickle.load(f) return peopleData, statsData def plot(people, stats): # experimental plotting function, currently need a rework using subplots numplots = len(stats.getStats()) colormap = plt.cm.gist_ncar plt.gca().set_color_cycle([colormap(i) for i in np.linspace(0, 0.9, numplots)]) labels =  for stat, label in stats.getStats(): plt.plot(stat) labels.append(label) plt.grid(True) plt.legend(labels, ncol=2, loc='upper center', bbox_to_anchor=[0.5, 1.1], columnspacing=1.0, labelspacing=0.0, handletextpad=0.0, handlelength=1.5, fancybox=True, shadow=True) # plt.yscale('log') plt.show() # Helpful diagnostic histogram for various stats # bins = np.arange(0, 100+1, 1) # plt.hist([person.maxAge for person in people], bins=bins) # plt.show() if __name__ == '__main__': # deserializing usually happens here, but I have excluded that portion here main()
I noticed that I needed a way to keep track of stats over time, so I could see this as plots and not just printed stats. I wrote a helper class
Stats to do this:
# stats.py class Stats(object): def __init__(self, oldDate, startingPopulation, statYear): self.timeTaken =  self.population =  self.meanAge =  self.newPeople =  self.deaths =  self.deltaNewPeople =  self.surnames =  self.deltaDeaths =  self.netPopDelta =  self.startingPopulation = startingPopulation self.year = 0 self.statYear = statYear self.oldDate = oldDate self.newPersonCount = 0 self.oldDeaths = 0 self.oldPeople = 0 def getStats(self): # convenience function to get all plot-able stats as an iterable with their labels attached return [ (self.timeTaken, "Time Taken"), (self.population, "Pop Count"), (self.meanAge, "Mean Age"), (self.newPeople, "New People"), (self.deaths, "Deaths"), (self.deltaNewPeople, "Delta New People"), (self.deltaDeaths, "Delta Deaths"), (self.surnames, "Surnames"), (self.netPopDelta, "Net Pop Change") ]
After a lot of trial and error and diagnostic plots, I feel I finally got this working decently. By no means do I claim that the results are statistically accurate for any real world population, but for what I set out to explore, the results seem reasonably OK, especially since The Atlantic notes that for China, with 1.3 billion people, 87% of the population shares only 100 unique surnames.
For a country of 1.3 billion people, there is a remarkably small number of common last names in China. An estimated 87 percent of the population shares one of 100 surnames, and more than one in five Chinese citizens is surnamed Li, Wang, or Zhang: more than 275 million people in all.
A sample run was done, simulating 500 years, printing stats every 100 years:
2018-03-08 15:59:02.250965 Time taken: 0 seconds (0 milliseconds / year) Current population size at year 0: 10000 Current mean age: 0 Current new people: 0 Current deaths: 0 Current delta new people: 0 (0 people / year) Current delta deaths: 0 (0 deaths / year) Current surname count: 10000 Current partner percentage: 0 Net population per year: 0 2018-03-08 15:59:11.135088 Time taken: 8 seconds (88 milliseconds / year) Current population size at year 100: 26898 Current mean age: 25 Current new people: 48750 Current deaths: 31852 Current delta new people: 48750 (487 people / year) Current delta deaths: 31852 (318 deaths / year) Current surname count: 2044 Current partner percentage: 51 Net population per year: 168 2018-03-08 15:59:24.385102 Time taken: 13 seconds (132 milliseconds / year) Current population size at year 200: 32645 Current mean age: 25 Current new people: 113218 Current deaths: 90573 Current delta new people: 64468 (644 people / year) Current delta deaths: 58721 (587 deaths / year) Current surname count: 1085 Current partner percentage: 52 Net population per year: 57 2018-03-08 15:59:41.181756 Time taken: 16 seconds (167 milliseconds / year) Current population size at year 300: 39817 Current mean age: 25 Current new people: 191768 Current deaths: 161951 Current delta new people: 78550 (785 people / year) Current delta deaths: 71378 (713 deaths / year) Current surname count: 806 Current partner percentage: 53 Net population per year: 71 2018-03-08 16:00:02.352846 Time taken: 21 seconds (211 milliseconds / year) Current population size at year 400: 47657 Current mean age: 25 Current new people: 286314 Current deaths: 248657 Current delta new people: 94546 (945 people / year) Current delta deaths: 86706 (867 deaths / year) Current surname count: 659 Current partner percentage: 50 Net population per year: 78 2018-03-08 16:00:28.200274 Time taken: 25 seconds (257 milliseconds / year) Current population size at year 500: 58446 Current mean age: 25 Current new people: 400610 Current deaths: 352164 Current delta new people: 114296 (1142 people / year) Current delta deaths: 103507 (1035 deaths / year) Current surname count: 578 Current partner percentage: 53 Net population per year: 107
Now to why I'm here on Code-Review in the first place with this monster. I feel I have whittled down some of the major performance issues, but fresh eyes tend to see things that I miss. I'm also looking to make things more Pythonic and readable. One thing I originally forgot to mention explicitly, though it does appear in the code, is that I want to have a more robust plotting scheme, probably utilizing subplots. I haven't quite gotten around to making that happen, so if you see a better way, please share!
I don't claim to be a python guru by any means, but StackOverflow goes a long way to help learn a new language. Feel free to discuss anything you find odd, and if anything is unclear as to what it's doing, feel free to ask.
In terms of readability, I am aware that one of the things
PEP8 specifies is
snake_case but I'm originally a c++ developer and old habits die hard (and I hate typing underscores). I also added a bunch of comments explaining my rationale behind why things appear as they do, which I feel gums things up a little, but I understand that everyone is not in my head.
Some side notes about timings below:
$ python -m timeit -s "import numpy as np" "np.random.randint(0, 1000)" 10000000 loops, best of 3: 0.164 usec per loop $ python -m timeit -s "import random" "random.randint(0, 1000)" 1000000 loops, best of 3: 1.52 usec per loop $ python -m timeit -s "import random" "random.uniform(0, 1)" 1000000 loops, best of 3: 0.489 usec per loop $ python -m timeit -s "import numpy as np" "np.random.uniform(0, 1)" 1000000 loops, best of 3: 0.219 usec per loop timeit.timeit("rand(0, 1000)", setup="from numpy.random import uniform as rand") 0.1250929832458496 timeit.timeit("rand(0, 1000)", setup="from random import uniform as rand") 0.47645998001098633