Demographic model

The original paper used a very simple instant population growth model. Under the model assumption, a population with an initial population size would evolve generations, instantly expand its population size to and evolve another generations. Such a model can be easily implemented as follows:

def ins_expansion(gen):
    'An instant population growth model'
    if gen < G0:
        return N0
    else:
        return N1

Other demographic models could be implemented similarly. For example, an exponential population growth model that expand the population size from to in generations could be defined as

def exp_expansion(gen):
    'An exponential population growth model'
    if gen < G0:
        return N0
    else:
        rate = (math.log(N1) - math.log(N0))/G1
        return int(N0 * math.exp((gen - G0) * rate))

That is to say, we first solve from and then calculate for a given generation.

There is a problem here: the above definitions treat N0, G0, N1 and G1 as global variables. This is OK for small scripts but is certainly not a good idea for larger scripts especially when different parameters will be used. A better way is to wrap these functions by another function that accept N0, G0, N1 and G1 as parameters. That is demonstrated in Example reichDemo where a function demo_model is defined to return either an instant or an exponential population growth demographic function.

Example: A demographic function producer

>>> import simuPOP as sim
>>> import math
>>> def demo_model(model, N0=1000, N1=100000, G0=500, G1=500):
...     '''Return a demographic function
...     model: linear or exponential
...     N0:   Initial sim.population size.
...     N1:   Ending sim.population size.
...     G0:   Length of burn-in stage.
...     G1:   Length of sim.population expansion stage.
...     '''
...     def ins_expansion(gen):
...         if gen < G0:
...             return N0
...         else:
...             return N1
...     rate = (math.log(N1) - math.log(N0))/G1
...     def exp_expansion(gen):
...         if gen < G0:
...             return N0
...         else:
...             return int(N0 * math.exp((gen - G0) * rate))
...     if model == 'instant':
...         return ins_expansion
...     elif model == 'exponential':
...         return exp_expansion
...
>>> # when needed, create a demographic function as follows
>>> demo_func = demo_model('exponential', 1000, 100000, 500, 500)
>>> # sim.population size at generation 700
>>> print(demo_func(700))
6309

now exiting runScriptInteractively...

Download reichDemo.py

Note

The defined demographic functions return the total population size (a number) at each generation beacuse no subpopulation is considered. A list of subpopulation sizes should be returned if there are more than one subpopulations.