gistfile1.py

# -*- coding: utf-8 -*-
import simuPOP as sim
import numpy as np
import random
from math import exp
#from scipy.stats import poisson # for poisson pmf

popHR = 2000
popSI = 1000
population_size = popHR + popSI

primiparity = 3 

def logismod(x, K, P0, k):
    res = K/(P0 + (K-P0) * exp(1 - k*x))
    return res

P0 = 0.1
K = 10
k = 0.7

pop = sim.Population(size = [popHR, popSI], loci = 3, infoFields = ["age", "ywc", "fitness", "hc", "iscub", "father_idx", "mother_idx", "ind_id"])
# ywc = years with cubs

random.seed(357)

# initialize infoFields, sex and genotype
sim.initInfo(pop = pop, values = np.random.poisson(3, population_size).tolist(), infoFields = "age")
sim.initInfo(pop = pop, values = 0, infoFields = "iscub")
sim.initInfo(pop = pop, values = 0, infoFields = "ywc")
sim.tagID(pop = pop)

sim.initSex(pop, maleProp = 0.5)

sim.initGenotype(pop = pop, prop = [0.1]*10)

fam_size = [[], []] 
fam_size_it = [[], []] 


def calculatePopulationSize(pop):
    global fam_size
    fam_size = [[], []]

    for sub in range(0, len(pop.subPopSizes())):
        # for each female sample litter size and store it in 'x'
        x = []
        for ind in pop.individuals(subPop = sub):
            if ind.sex() == sim.FEMALE and ind.hc == 0 and ind.age >= primiparity:
                x.append(np.random.choice(a = [1, 2, 3, 4], size = 1, p = [0.15, 0.65, 0.15, 0.05], replace = True).tolist())
        fam_size[sub] = sum(x, [])

    # sum population size and extra cubs population-wise
    post_rep_popsize = []
    for i in range(0, len(pop.subPopSizes())):
        # print("[calculatePopulationSize] subPopSizes: {} added {} cubs (total {})".format(pop.subPopSizes()[i], sum(fam_size[i]), pop.subPopSizes()[i] + sum(fam_size[i])))
        post_rep_popsize.append(pop.subPopSizes()[i] + sum(fam_size[i]))

    # make 'fam_size' a flat list
    fam_size = sum(fam_size, [])
    # print("calculatePopulationSize: {}".format(fam_size))

    global fam_size_it
    fam_size_it = fam_size
    # return final population size which consists of cloned non-offspring
    # and new offspring
    #print("calculatePopulationSize: population size {}".format(post_rep_popsize))
    # print("===")
    return post_rep_popsize

# pick random father
def bearFather(pop, subPop):
    # define fitness for individuals based on a logistic model
    # see parametrization at the beginning of the file
    all_males = [x for x in pop.individuals() if x.sex() == sim.MALE]
    for ind in all_males:
        ind.fitness = logismod(x = ind.age, K = K, P0 = P0, k = k)

    # prepare males
    allMales = sim.RandomParentChooser(sexChoice = sim.MALE_ONLY,
                                       selectionField = "fitness",
                                       replacement = True)
    allMales.initialize(pop, subPop) # builds individual index

    while True:
        male = allMales.chooseParents()[0]
        # print(male.fitness)
        yield male

# pick a mother that has no cubs
def bearMother(pop, subPop):
    while True:
        # prepare females in estrus
        all_females = [x for x in pop.individuals(subPop) if x.sex() == sim.FEMALE and x.hc == 0 and x.age >= 3]
        #print("bearMother: number of females is {}".format(len(all_females)))
        pick_female = np.random.randint(0, (len(all_females)))
        female = all_females[pick_female]
        female.hc = 1
        female.ywc = 0
        
        yield female

# Read family size from a global variable and use that as
# the number of offspring per mating event. This function is applied
# to individual subpopulation
def familySizeGenerator():
    global fam_size_it

    # print("familySizeGenerator: === population ===")
    # print("familySizeGenerator: fam_size_it = {}".format(fam_size_it))
    # print("familySizeGenerator:    fam_size = {}".format(fam_size))
    for sz in fam_size_it:
        # print("familySizeGenerator: producing {} cubs".format(sz,))
        fam_size_it = fam_size_it[1:] # after cub has been "made", remove it
        yield sz

# debugging functions
def endGenerationClick():
    print("=== END OF GENERATION TALLY ===")
    print("fam_size: {}". format(fam_size))
    print("fam_size_it: {}". format(fam_size_it))
    print("=== END OF GENERATION TALLY ===")
    return True

def updateYWC(ind):
    # cubs become independent of their mother after 'primiparity'
    if ind.hc == 1 and ind.ywc >= primiparity:
        ind.hc = 0
        ind.ywc = 0

    # if female is rearing cubs, add 1 to 'ywc' (years with cubs)
    if ind.hc == 1:
        # print("[updateYWC] YWC: {}". format(ind.ywc))
        ind.ywc += 1

    # if individual is '>= primiparity' release it from its mother
    if ind.iscub == 1 and ind.age >= primiparity:
        ind.iscub == 0

    return True

def cullSlovenia(pop, param):
    # param is a dictionary {key: value}
    # lam (numeric) Lambda parameter in Poisson process, e.g. {"lam": 1.82}

    # generate expected values according to the distribution
    raw_cull = np.random.poisson(lam = param["lam"], size = param["cull_size"]).tolist()
    hist = [[x, raw_cull.count(x)] for x in set(raw_cull)]

    # "tabulate" values into a histogram, beware, object is a dictionary
    tabulated_cull = {x[0]: x[1] for x in hist}

    print("CULLSLOVENIA: calculated frequencies by age classes:")
    print(tabulated_cull)

    # find all individuals in subpopulation of the said age class
    for n, v in tabulated_cull.iteritems(): # n = class, v = frequency
        n_before = pop.popSize()
        inds = [x for x in pop.individuals() if x.age == n and x.hc == 0]

        # sample out v number of individuals
        inds = random.sample([x for x in inds], v)

        # use ind_ids of individuals to remove them from the (sub)population
        pop.removeIndividuals(IDs = [x.ind_id for x in inds], idField = "ind_id")
        n_after = pop.popSize()

        print("CLASS [{}]: BEFORE [{}] AFTER [{}] CHANGE [{}]".format(n, n_before, n_after, (n_before - n_after)))

    return True

pop.evolve(
    preOps = [
        # sim.IdTagger(infoFields = "ind_id"),
        sim.Stat(popSize = True),
        sim.InfoExec('age += 1'), # increase age for everyone
        sim.PyOperator(func = updateYWC) # bookkeeping for mothers with cubs
    ],
    matingScheme = sim.HeteroMating([
                                        sim.HomoMating(
                                            chooser = sim.CombinedParentsChooser(
                                                fatherChooser = sim.PyParentsChooser(generator = bearFather),
                                                motherChooser = sim.PyParentsChooser(generator = bearMother)
                                            ),
                                            generator = sim.OffspringGenerator(ops = [
                                                sim.InfoExec("age = 0"),
                                                sim.InfoExec("iscub = 1"),
                                                sim.ParentsTagger(),
                                                sim.IdTagger(),
                                                sim.MendelianGenoTransmitter()
                                            ] ,
                                                                               numOffspring = familySizeGenerator)),
                                        sim.CloneMating(subPops = sim.ALL_AVAIL, weight = -1)
                                    ],
                                    shuffleOffspring = False,
                                    subPopSize = calculatePopulationSize
    ),
    postOps = [
        #this line has been left empty on purpose
        sim.PyOperator(func = cullSlovenia(pop, param = {"lam": 1.82, "cull_size": 50}), subPops = 1),
    ],
    gen = 3
)