BoPeng/simuPOP

Handle cases where no males/females exist in a *virtual* subpopulation?

ashander opened this issue · 35 comments

I'm attempting to deal with virtual subpopulations where there are no male or females using an approach similar to your suggestion in #15 and running into problems (related to weighting; see also #12 ).

Essentially, I can control the next generation size (as here with demoModel) and even the offspring number and sequencing (although here that is static) but when multiple virtual subpopulations exist in the same physical subpop, I'm running into issues due to weighting --- the virtual subpopulation with no females still attempts to produce offspring because weights are assigned based on current size (see debug output below) and so assigns one offspring to each VSP (VSP sizes in subpop 0 is 1, 1).

Is there some weighting scheme (or other strategy) that will work to assign zero offspring to VSPs with no males/females?

import simuPOP as sim
numoff = 1  # Global used in in demoModel and evolve

def demoModel(pop):
    sim.stat(pop, numOfMales=True, vars=['numOfMales_sp', 'numOfFemales_sp'])
    m_f = [(pop.dvars(i).numOfMales, pop.dvars(i).numOfFemales) for i in range(pop.numSubPop())]
    return [f * numoff if m > 0 else 0 for m,f in m_f]

sim.turnOnDebug('DBG_ALL')
sim.setOptions(seed=112)
initial_size = 7
p = sim.Population(initial_size, infoFields=['mark'])
sim.initInfo(p, [1, 1, 1, 2, 2, 2, 2], infoFields='mark')
sim.initSex(p, sex=[sim.MALE, sim.FEMALE, sim.FEMALE, sim.MALE, sim.MALE, sim.MALE, sim.MALE, ])
p.setVirtualSplitter(sim.InfoSplitter('mark', ranges=((1, 2), (2, 3))))
p.evolve(
    preOps=sim.Dumper(structure = False),
    matingScheme=sim.HeteroMating(
        [sim.HomoMating(
            chooser = sim.CombinedParentsChooser(
                sim.RandomParentChooser(replacement=True, sexChoice=sim.MALE_ONLY),
                sim.SequentialParentChooser(sexChoice=sim.FEMALE_ONLY),
            ),
            generator = sim.OffspringGenerator( ops = [
                sim.MendelianGenoTransmitter(),
            ],
                numOffspring=numoff),
            subPops=[(sim.ALL_AVAIL, 0), (sim.ALL_AVAIL, 1)],
            weight=1,
        )],
        subPopSize=demoModel,
    ),
    postOps=sim.Dumper(structure = False),
    gen=1,
)

Output

::::::::::::::::::::::::::::::::--
Constructor of population is called
Inc ref of 2 to 9
individual size is 24+1*0
, infoPtr: 8, GenoPtr: 8, Flag: 1, plus genoStru
genoSize 0
Analyzed string is >
Filename is 
Analyzed string is >
Filename is 
Creating Simulator 
Constructor of population is called
Inc ref of 1 to 6
individual size is 24+1*0
, infoPtr: 8, GenoPtr: 8, Flag: 1, plus genoStru
genoSize 0
Constructor of population is called
Inc ref of 1 to 7
individual size is 24+1*0
, infoPtr: 8, GenoPtr: 8, Flag: 1, plus genoStru
genoSize 0
Simulator created
Start evolution.: 2e-06
SubPopulation 0 (), 7 Individuals:
   0: MU  |  |  1
   1: FU  |  |  1
   2: FU  |  |  1
   3: MU  |  |  2
   4: MU  |  |  2
   5: MU  |  |  2
   6: MU  |  |  2

Applied <simuPOP.Dumper>: 0.000186
Start mating at generation 0: 1.3e-05
Dec ref of 1 to 6
Inc ref of 2 to 10
New subpop size 2
Positive mating scheme weights: 1, 1
Negative mating scheme weights: 0, 0
VSP sizes in subpop 0 is 1, 1
Mating is done in single-thread mode
Mating scheme 1 in subpopulation 0 failed to produce 1 offspring.

It seems that I have somehow missed this ticket. Let me know if you still need help. Basically, the weight design does not consider proportion of males and females within VSPs so the only way to avoid this is to check VSPs subpopulation before mating.

Thanks! I should've replied that we solved it another way (with a custom parent chooser that does essentially what you're suggesting).

But, it does seem like a somewhat fundamental issue that many folks may run into. Are there large obstacles to a low-level patch that would correct the weight design for cases like my example above? If you don't see any, I could try to contribute one =)

There are a few problems that make it difficult to produce a general solution:

  1. Even if VSP does not have any male or female, mating can still happen if a non-sexual mating scheme is used, or if a conditional mating scheme is used to perform selfing in such cases to mimic the behavior of some plant populations. Basically, such cases can be intentional.
  2. We cannot delay the problem to the mating stage because a pre-specified number of offspring (calculated from weight) has to be generated to populate the offspring subpopulation.

So an option is needed here and it has to determine how the size of offspring populations is determined. For example,

  1. weight_by=POPSIZE: current behavior
  2. BY_MALE: only count male individuals in the parental subpopulation
  3. BY_FEMALE: only count female individuals in the parental subpopulation
  4. BY_PAIR: determine the number of "mating" parents by number of pairs, namely min(male, female).

The last option would treat SP (or VSP) with no male or no female as having zero mating parents and produce zero offspring. However, it can be confusing to use BY_MALE with weight=-1 because weight=-1 was used to generate the same number of offspring as parents, and now it would be the same number of offspring as fathers (so weight=-2 makes more sense).

OK, this is what has been implemented:

As a special case that can be quite annoying during the simulation of small populations, a (virtual) subpopulation can have no male and/or female. If the parental (virtual) subpopulation is empty, it will produce no offspring regardless of its weight. However, if the parental (virtual) subpopulation is not empty, it will be expected to produce some offspring, which is not possible if a sexual mating scheme is used. In this case, you can use a parameter weightBy to specify how parental (virtual) population sizes are calculated. This parameter accepts values ANY_SEX (default), MALE_ONLY, FEMALE_ONLY, PAIR_ONLY, and use all individuals, number of male individuals, number of female individuals, and number of male/female pairs (basically the less of numbers of males and females) as the size of parental (virtual) subpopulation, respectively. When weightBy=PAIR_ONLY is used, parental (virtual) subpopulations with only males or females will appear to be empty and produce no offspring. Note that in this mode (also MALE_ONLY and FEMALE_ONLY), the perceived parental population sizes are no longer the actual parental population sizes so you might need to adjust parameter weight (e.g. weight=-2) to produce correct number of offspring.

Hello,
I feel like the parameter weightBy is intended to be optional, though I cannot manage to make my code work without including it (i.e. as it was with the 1.1.8.3 version of simuPOP). Is it meant to be mandatory ?
Thank you for your answer.

Yes, weightBy is meant to be optional and has a default value ANYSEX (f5a22bb#diff-b43b9f8d4f6b1e382f8e2a7c160c1921R1730). Could you provide an example where weightBy is required? Note that the latest version of simuPOP is 1.1.9 but this should not affect this particular feature.

At least the basic HeteroMating example https://github.com/BoPeng/simuPOP/blob/master/docs/HaplodiploidMating.py works without this option.

Indeed, this example works also on my setup.
This is how the mating scheme is defined in my code :

mating = simu.HeteroMating([
            simu.SelfMating(ops=[
                simu.Recombinator(rates=self.metaPop.calcRecombinaison(), 
                             loci=range(len(self.metaPop.getMarkers())-1))],
                 weight=self.metaPop.getSelfing()),
            simu.HermaphroditicMating(ops=[
                   simu.Recombinator(rates=self.metaPop.calcRecombinaison(), 
                       loci=range(len(self.metaPop.getMarkers())-1))], 
                weight=1-self.metaPop.getSelfing())
              ], 
         subPopSize=self.reproduction)

I don't know how to adapt it to this example... And not sure this is actually relevant to this thread.
I hope everything is quite clear with this code. Let me know if I should further describe it.
The error I get (and which I had not with version 1.1.8.3) is : ValueError: mating.cpp: 1703 Overall weight is zero

Could you run

turnOnDebug('DBG_DEVEL')

before the evolve function and report the output of this debug message?

So, did you say the problem would disappear if you added weightBy somewhere to your example? How did you add it?

The output of the debug function before the error is :
VSP sizes in subpop 55 is 0, 10 Parental Population Size: 0, 0 Positive mating scheme weights: 0, 0 Negative mating scheme weights: 0, 0

Sorry if I wasn't clear with what I meant... I was previously using v1.1.8.3 which had not the weightBy parameter and was working. In the new version (1.1.9 on my system) my mating doesn't work (as shown above). I was assuming the problem was related to this parameter. I'm still working on using this parameter in my code though I might need quite a lot of time.

The debug message says that the first parental population is empty when the problem happen. I can reproduce your problem with the following example

import simuPOP as simu

pop = simu.Population(size=[0, 55], loci=2)

mating = simu.HeteroMating([
    simu.SelfMating(ops=[
                simu.Recombinator(rates=0.001, loci=0)],
        weight=0.5),
    simu.HermaphroditicMating(ops=[
               simu.Recombinator(rates=0.1, loci=0)], 
        weight=0.5)
    ])

pop.evolve(
        initOps=simu.InitSex(),
        matingScheme=mating,
        gen=1)

but I do not see how weightBy would help...

Sorry, I might have been mislead by the changelog notifying this change in the function I was using...
Any idea why 1.1.8.3 would deal with size 0 populations but not 1.1.9 ?

This can be a regression bug but you seemed to say that you can use weightBy to fix the problem. If this is the case, could you let me know how to fix the problem by adding weightBy? If not, you will have to revert to 1.1.8.3 before I fix it and release 1.1.9.1.

I tried several things to make it work with weightBy, thought it didn't work. I'll revert to v.1.1.8.3 until v1.1.9.1 is released. Thank you for your help. Do I need to do anything (create an issue or whatever) ?

Thanks. This is good enough. I will reopen the ticket and work on it hopefully later today.

Hi,

I have run into the same problem while running captive breeding simulations (RuntimeError: RandomParentsChooser fails because there is no female individual in a subpopulation). The species I am simulating is a protandric hermaphrodite (starts male and switches to female, and can also switch back and forth throughout its life), so a function that calculates males and females and changes the sex accordingly before or during mating would be ideal. However, I am a python novice and am having trouble coding this function. I know where to put it in the sims, but was wondering if you could share how you coded a sex change? Thank you in advance.

Cheers,
Katie

Sex change is easy, something like

if num_male == 0: # need to change some of the parents to female
    for i in range(num_to_change):
        pop.individual(i, subpop_index_or_ignore_for_all).setSex(sim.FEMALE)

This piece of code can be used in a demographic function or pre-mating PyOperator so that it can be called right before mating happens.

This is a huge help and I really appreciate it. Your prompt replies and willingness to help along with the detailed simuPOP manual have made it possible for a python novice like myself to make sufficient progress with coding complex simulations, so thank you so much for that.

Hi,
It seemed as though the sex change function was working, but I just decided to test it to be sure by setting the population to all females and now it is not working (so I guess it never really was working). It is probably a simple fix that I am missing–it seems as though the stat is not passing to numOfMales.
Here is my function:

def sexm(pop):
    for i in range(50):
        pop.individual(i).setSex(sim.MALE)
        return True

Here is when I call it in the preOps:

sim.Stat(pop, numOfMales = True),
        sim.IfElse('numOfMales==0',
               ifOps = sim.PyOperator(func = sexm)
               ),

error:

RuntimeError: RandomParentsChooser fails because there is no male individual in a subpopulation.

If you are not sure if the IfElse stuff works, you can bypass it by calling the following function in PyOperator.

def sexm(pop):
    sim.stat(pop, numOfMales=True)
    print(f'The number of males is {pop.dvars().numOfMales}')
    if pop.dvars().numOfMales == 0:
        for i in range(50):
            pop.individual(i).setSex(sim.MALE)
    return True

Thank you for your prompt responses and help. The strangest thing is that both functions (yours above and mine below to include vsps since I have age structure and only ages 2-10 mate) fail when I initialize with either 0 males or 0 females (I have functions for both change from male to female and vice versa). I'm thinking this may be a result of my mating scheme (maybe weight?) or the function isn't working at all? For example, I initialized with a small # of individuals (100) and had it print out the num of males and females in each vsp:

wild - age < 1: Females 0 Males 0 
wild - 1 <= age < 2: Females 2 Males 5 
wild - 2 <= age < 3: Females 0 Males 2 
wild - 3 <= age < 4: Females 3 Males 1 
wild - 4 <= age < 5: Females 2 Males 4 
wild - 5 <= age < 6: Females 1 Males 9 
wild - 6 <= age < 7: Females 2 Males 5 
wild - 7 <= age < 8: Females 4 Males 5 
wild - 8 <= age < 9: Females 1 Males 7 
wild - 9 <= age < 10: Females 1 Males 1 
wild - 10 <= age < 11: Females 2 Males 5 

Overall:

The number of males is 39
The number of females is 16

Error:

Mating scheme 12 in subpopulation 0 failed to produce 19 offspring.
Traceback (most recent call last):
RuntimeError: RandomParentsChooser fails because there is no female individual in a subpopulation 

I get the same error when I initialize with no males as well and it prints:
The number of males is 0 and the same error follows.

Here is my function for females:

def sexf(pop):
    sim.stat(pop, numOfMales=True)
    print(f'The number of females is {pop.dvars().numOfFemales}')
    if pop.dvars().numOfFemales == 0:
        for i in range(pop.numVirtualSubPop()):
            pop.individual(i).setSex(sim.FEMALE)
    return True

Here is my mating scheme:

def demo(pop):
    return [x + random.uniform(200, 600) for x in pop.subPopSizes()]

matingScheme=sim.HeteroMating([
        # CloneMating will keep individual sex and all
        # information fields (by default).
        sim.CloneMating(subPops=[(0,sim.ALL_AVAIL)], weight=-1),
        # only individuals with age between 2 and 10 will mate
        #  The age of offspring will be zero.
        sim.RandomMating(ops=[
            sim.IdTagger(),  # give new born an ID
            sim.PedigreeTagger(),  # track parents of each individual
            sim.MendelianGenoTransmitter(),  # transmit genotype
        ],
            numOffspring=(sim.GEOMETRIC_DISTRIBUTION, 0.01),
            subPops=[(0,2),(0,3),(0,4),(0,5),(0,6),(0,7),(0,8),(0,9),(0,10)]),],
        subPopSize=demo),

I somehow missed this ticket. @khornick1 , do you still have trouble with it?

No problem. To deal with this issue, I just forced the reproducing vsps (ages 2-11) to ~50/50 male and female, which actually describes the captive breeding protocols of the oyster hatchery I am trying to simulate (and the species is a protandric hermaphodrite so it changes sex). Thanks for checking, though I appreciate it.

Actually I'm getting this error again, did you have any advice? Thanks.

From your original post, you have zero male and female in vsp0, and that is why you got the error. What did you do to get around of the problem? Could you post a complete script that demonstrate the problem? Please also upgrade to the latest version of simuPOP (through conda) if you are not using it so that we are on the same page.

My simuPOP version is all up to date. That's strange though, because vsp0 doesn't mate so it shouldn't matter whether there are males and females for the randomparentchooser? Before my workaround was to "force" sexes prior to mating scheme using sim.InitSex. Now that I have changed my mating scheme around a bit and am using different numbers, I am met with the same error and haven't had luck with implementing a sex change function to get around it. The error I get (when running through jupyter notebooks):

Mating scheme 21 in subpopulation 2 failed to produce 8 offspring.
RuntimeError Traceback (most recent call last)
in
263 sim.PyOutput('\n', step=1),
264 ],
--> 265 gen=pre_hatchery_gens
266 )
267

~/miniconda3/envs/simupop/lib/python3.6/site-packages/simuPOP/init.py in evolve_pop(self, initOps, preOps, matingScheme, postOps, finalOps, gen, dryrun)
405 simu = Simulator(self)
406 # evolve
--> 407 gen = simu.evolve(initOps, preOps, matingScheme, postOps, finalOps, gen)
408 # get the evolved population
409 self.swap(simu.population(0))

RuntimeError: RandomParentsChooser fails because there is no female individual in a subpopulation

I have attached my script which only includes the first of three phases of my simulations, of which the first I am getting the error above. Thank you in advance for your help.
simsaug2019.txt

It took a while but I could reproduce the problem. I will have a closer look later today.

Hi Bo,

Thanks for your time and help. I think the error occurs because one of vsps does not have a male or female in it. Perhaps some sort of check prior to mating to make sure each vsp has a male or female in it might solve the issue? I guess I thought if there were males in females in the mating population, then mating should happen but this suggests that mating occurs only within individuals in the same vsp?

As far as I can tell, the failure happens at the third mating scheme (RandomMating) that is applied to VSPs of different ages, and then one of the age group has no male or female.

This thread was created because of similar problems and a solution was to use option weightBy, which in your case would be "weightBy=PAIR_ONLY", meaning if a VSP has only male or female, it would have weight zero.

I added weightBy before option subPopSize as follows

        weightBy=sim.PAIR_ONLY,
        subPopSize=demo),

and the simulation has been running for a while. Note that weight by pair will lead to less than half parental population size (e.g. a VSP with 8 males and 5 females will have a parental population size of 5), so you should at least change the weight of CloneMating from -1 to -2, and not everyone will be copied if you do not have matching number of males and females (e.g. 10 offspring will be copied in the VSP).

Okay, got it. Thank you so much for your help and for taking the time to run my script.

Just wondering, could another solution be to check any of the mating VSPs and if one of them lacks a male or female change the sex prior to mating?

Of course, the demographic function will be called right before mating and is the right place for this kind of thing. Basically you can apply sim.Stat(num_of_males) to the VSPs, get the number of males and females and change sex if needed.

Ugh, I just realized a bit late that I would like mating to occur between different age classes, for example, a 2 year old can mate with an 11 year old (which also might help with this error). What would be the easiest way to do this without changing too much in my code? Define another type of splitter for mating individuals?

You can adjust your VSP to be "coarser", not by each year, but by range of years, something "junior", "adult", and "senior", and adult would be from 2 to 11. If you defined VSPs for each year, you can add new VSPs to existing ones (see
CombinedSplitter in http://simupop.sourceforge.net/manual_svn/build/refManual_ch2_sec2.html#class-combinedsplitter). When VSPs get complex, you can use names instead of numbers to avoid mistake (e.g. [0, "adult']).

Thank you!