BoPeng / simuPOP

A general-purpose forward-time population genetics simulation environment.
http://bopeng.github.io/simuPOP/
GNU General Public License v2.0
31 stars 12 forks source link

Variable Age at First Breeding #36

Closed snunzi closed 7 years ago

snunzi commented 7 years ago

Hi Everyone, I am attempting to build a model where individuals reach breeding age at 1 to 5 years old. I am trying to set it up so that x% return for the first time in Year 1, y% in Year 2, etc. up through Year 5, and the sum of those first five years of first-time breeders equals 100%. Once an animal does return for the first time, then it has a set probability of breeding from then on. Do you think simuPOP is capable of setting up a model like this? Would I have to have a VSP for each cohort born? It seems like this would get very complicated with overlapping generations, maybe there is an easier way I am overlooking? Any insights on this would be much appreciated!

Thanks, Schyler

Breeding_Probability.xlsx

BoPeng commented 7 years ago

This can be done through VSPs but in a complicate scenario like this, I would use a PyParentsChooser to explicitly choose parents for breeding. Basically you create groups of parents according to age and choose parents by certain probabilities. The code is longer but you have control on exactly who would be chosen. Note that you can set some flag for mother who have produce offspring (perhaps an during-mating PyOperator) if a mother cannot breed again within the first five years.

snunzi commented 7 years ago

Thank you, Bo. Yes, it makes much more sense to do it like that!

One more thing that I'm stuck on is how to assign a cohort density to an offspring. I have an infofield called 'cohort_density' and basically I want to tag each offspring as coming from either a low or high density cohort. So, if the total number of newborns is > 10,000 it would be off.setInfo(0, 'cohort_density'), and if the total number of newborns is < 10,000 it would be off.setInfo(1, 'cohort_density'). I think I should be able to do a simple PyOperator during mating like below, but can't find a code that works. Thank you so much again!

mating_scheme = sp.HeteroMating([
    sp.HomoMating(
        sp.CombinedParentsChooser(
            fatherChooser = sp.PyParentsChooser(FatherChooser) ,
            motherChooser = sp.PyParentsChooser(MotherChooser)),
        sp.OffspringGenerator(numOffspring=famSizeGenerator, ops=[
            sp.InfoExec("age = 0"),
            sp.MendelianGenoTransmitter(),
            sp.IdTagger(),
            sp.ParentsTagger(),
            sp.PyOperator(lambda mom: mom.setInfo(1, 'has_offspring') is None),
            sp.PyOperator(lambda dad: dad.setInfo(1, 'has_offspring') is None), 
            sp.PyOperator('cohort_density' ????)             ]),
        weight=1),
    sp.CloneMating(subPops = sp.ALL_AVAIL, weight=-1),
    ],
    subPopSize=demoModel)
BoPeng commented 7 years ago

I do not quite get it. Is cohort_density determined by the number of newborns in the population? If that is the case, this would be a population level info, not individual level, and logically speaking cannot be determined during mating, right?

Or, is this something you would like to derive (inherit) from parents? In that case you would assign this info to parents using some criteria and then some operator like Inherittagger to pass from parent to offspring.

snunzi commented 7 years ago

I think I have figured it out. I determine the total number of newborns added to the population by calling demoModel function during mating. Since the number of newborns is already determined, this information is used in the cohortdensity function to assign the newborns the cohort_density infoField (based on the total cohort size in which they were born). I am wondering if you could skim my cohortdensity function below and see if you think it makes sense. Thank you again.

-Schyler

#This function chooses parents for mating according to age-based fecundity
rates, and limits number of offspring mothers can produce per breeding
cycle according to their age
famSize = []
mothers = []
fathers = []

def demoModel(pop):
    global famSize
    global mothers
    global fathers
    famSize = []
    mothers = []
    fathers = []
    for ind in pop.individuals():
        if ind.sex() == 2 and ind.has_offspring == 0 and ind.cohort_density
== 1:
            mothers.extend([ind] * (random.random() <
pop.dvars().female_age_fecundity[int(ind.age)] ) )
        if ind.sex() == 2 and ind.has_offspring == 0 and ind.cohort_density
== 0:
            mothers.extend([ind] * (random.random() <
pop.dvars().female_age_fecundity_lowdens[int(ind.age)] ) )
        if ind.sex() == 2 and ind.has_offspring == 1:
            mothers.extend([ind] * (random.random() < 0.66 ) )
    for ind in pop.individuals():
        if ind.sex() == 1 and ind.has_offspring == 0 and ind.cohort_density
== 1:
            fathers.extend([ind] * (random.random() <
pop.dvars().male_age_fecundity[int(ind.age)] ) )
        if ind.sex() == 1 and ind.has_offspring == 0 and ind.cohort_density
== 0:
            fathers.extend([ind] * (random.random() <
pop.dvars().male_age_fecundity_lowdens[int(ind.age)] ) )
        if ind.sex == 1 and ind.has_offspring == 1:
            fathers.extend([ind] * (random.random() < 0.9 ) )
    for ind in mothers:
        if ind.cohort_density == 0:
            famSize.append(pop.dvars().max_kids_lowdens[int(ind.age)])
        else:
            famSize.append(pop.dvars().max_kids[int(ind.age)])
    meanoffspring=np.mean(famSize)
    print('Number of mothers: {} with {} eggs and mean per mom
{}'.format(len(mothers), sum(famSize), meanoffspring))
    print('Number of fathers: {}'.format(len(fathers)))
    if len(fathers) >= 1 and len(mothers) >= 1:
        return pop.popSize() + sum(famSize)
    else:
        return pop.popSize()

def famSizeGenerator():
    global famSize
    for sz in famSize:
        #print('producing {} eggs'.format(sz))
        yield sz

def MotherChooser():
    global mothers
    for ind in mothers:
        yield ind

def FatherChooser():
    global fathers
    while True:
        father = random.choice(fathers)
        yield father

def cohortdensity():
    global famSize
    if sum(famSize) >= 18700:
        sp.PyOperator(lambda off: off.setInfo(1, 'cohort_density') is None)
    else:
        sp.PyOperator(lambda off: off.setInfo(0, 'cohort_density') is None)
    return True

#Clone mating will move all parents on to the next generation, while
homomating will produce offspring according to parents chosen with father
and mother chooser
mating_scheme = sp.HeteroMating([
    sp.HomoMating(
        sp.CombinedParentsChooser(
            fatherChooser = sp.PyParentsChooser(FatherChooser) ,
            motherChooser = sp.PyParentsChooser(MotherChooser)),
        sp.OffspringGenerator(numOffspring=famSizeGenerator, ops=[
            sp.InfoExec("age = 0"),
            sp.MendelianGenoTransmitter(),
            sp.IdTagger(),
            sp.ParentsTagger(),
            sp.PyOperator(lambda mom: mom.setInfo(1, 'has_offspring') is
None),
            sp.PyOperator(lambda dad: dad.setInfo(1, 'has_offspring') is
None),
            sp.PyOperator(cohortdensity)]),
        weight=1),
    sp.CloneMating(subPops = sp.ALL_AVAIL, weight=-1),
    ],
    subPopSize=demoModel)
BoPeng commented 7 years ago

Your code looks complicated.... what does cohortdensity do? I mean,

def cohortdensity():
    global famSize
    if sum(famSize) >= 18700:
        sp.PyOperator(lambda off: off.setInfo(1, 'cohort_density') is None)
    else:
        sp.PyOperator(lambda off: off.setInfo(0, 'cohort_density') is None)
    return True

this function creates PyOperator but is not returned or applied. Are you trying to do something like the following?

def cohortdensity(off):
    global famSize
    off.cohort_density = sum(famSize) >= 18700
    return True
snunzi commented 7 years ago

Yes, it has grown very complicated. But your code help is so useful and appreciated. Ideally, cohortdensity function counts number of newborns and assigns each newborn cohort_density=1 if more than 18700 newborns are produced, and cohort_density=0 if less than 18700 newborns are produced. Given your code advice it would look like:

def cohortdensity(off): global famSize if sum(famSize) >= 18700: off.cohort_density = 1 else: off.cohort_density = 0 return True

On Thu, Apr 6, 2017 at 9:52 AM, Bo notifications@github.com wrote:

Your code looks complicated.... what does cohortdensity do? I mean,

def cohortdensity(): global famSize if sum(famSize) >= 18700: sp.PyOperator(lambda off: off.setInfo(1, 'cohort_density') is None) else: sp.PyOperator(lambda off: off.setInfo(0, 'cohort_density') is None) return True

this function creates PyOperator but is not returned or applied. Are you trying to do something like the following?

def cohortdensity(off): global famSize off.cohort_density = sum(famSize) >= 18700 return True

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/BoPeng/simuPOP/issues/36#issuecomment-292180906, or mute the thread https://github.com/notifications/unsubscribe-auth/AP53hruYVdjQDm499Qih7-kQI-vjHtfFks5rtO4LgaJpZM4M0UCM .

BoPeng commented 7 years ago

Isn't famSize determined before mating? If so, all offspring would have the same cohort_density, right? It then would be easier and perhaps more efficient to define a VSP for offspring (by age==0 I guess), and set this info globally in postOps, something like PyOperator(lambda pop: pop.setIndInfo(sum(famSize) > 18700, 'cohort_density', (0, 0)) is None) (if (0,0) is the VSP for newborn).

snunzi commented 7 years ago

That is more straight forward and works well. Thank you for your help.

On Thu, Apr 6, 2017 at 10:53 AM, Bo notifications@github.com wrote:

Isn't famSize determined before mating? If so, all offspring would have the same cohort_density, right? It then would be easier and perhaps more efficient to define a VSP for offspring (by age==0 I guess), and set this info globally in postOps, something like PyOperator(lambda pop: pop.setIndInfo(sum(famSize) > 18700, 'cohort_density', (0, 0)) is None) (if (0,0) is the VSP for newborn).

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/BoPeng/simuPOP/issues/36#issuecomment-292199665, or mute the thread https://github.com/notifications/unsubscribe-auth/AP53hnQQWA36x-PR06rxuIHkSd1pQCn2ks5rtPxigaJpZM4M0UCM .