exponential growth rates not working in fsc2702

[J Weir]

Dear Dr. Excoffier,

I cannot get your example (3PopExpBot20Mb) or my own models with an exponential growth parameter to work in fsc2702 (Linux), even though they work fine in fsc26. I suspect exponential growth has a bug in version 7. For example, I ran with the command line from the manual on the 3PopExpBot20Mb dataset as follows:

./fsc2702 -t 3PopExpBot20Mb.tpl -n100000 -d -e 3PopExpBot20Mb.est -M -L40 -q --multiSFS -C10 -c8

Some runs work, others crash with the following error:

Estimating model parameters using 12 batches and 8 threads
free(): invalid next size (fast)
Aborted (core dumped)

Given I cannot get your example with exponential growth to consistently work without crashing on at least some runs I suspect there could be a bug in the implementation of exponential growth in version 7.

As an aside, in all of the models I tested, I added the appropriate events to the TPL file to indicate a deme has zero growth rate (these were already included in TPL file for 3PopExpBot20Mb , though not shown in the manual itself), following its merger into another lineage. (e.g. "TDIV 1 1 0 1 0 1 "). I believe this event is a new requirement of fastsimcoal2 version 7, though its requirement is puzzling given I would have assumed that a line such as "TDIV 1 0 1 1 0 1" should be enough to tell fastSIMCOAL2 that deme 1 not only has a growth rate of 0, but no longer exists (since all of its population has migrated to deme 0  leaving deme 1 extinct). If it no longer exists then events like "TDIV 1 1 0 0 0 1" to kill the deme also seem unnecessary, but the manual seems to be promoting the use that line in certain contexts. I  find these new additions to the historical events difficult to understand.

Best wishes,

Jason Weir

[zhao wei]

I have a similar problem when using the Fsc27 to run exponential growth scenarios, The same scenarios work fine in fsc26. 

Best regards,

Wei

[Ivan Scotti]

Dear Jason and Wei,

in our recent experiences with Fsc27 on linux and MacOs machines, the example you used would work, but our own models would only work if the populations supposed to have "disappeared" were actively killed by adding a historical event, as indicated at page 29 of the manual. When we did not do it, the merged populations stayed around and continued to grow or shrink (we dubbed them zombie populations), eventually reaching infinity or zero size, thus making the program crash.

I hope this helps,

Regards

Ivan

[Laurent Excoffier]

Hi,

I think Scotti is right. Population size going to zero due to negative exponential growth causes problems. It is a feature, not a bug...

Note that populations always exist, even after all its lineages have been transferred to another deme. this is because you could have some migrants coming back to it later on (looking backward in time).

best

laurent

[J Weir]

Hi Laurent,

"Population size going to zero due to negative exponential growth causes problems. It is a feature, not a bug..."  Sure, but my understanding is that the way 3PopExpBot20Mb.tpl is coded it was supposed to have stopped negative exponential growth in demes 1 and 2 at TDIV before demes ever reached a size of 0. This is because R1 is coded such that there must always be at least 10 and no more than 10000 individuals in these demes at TDIV. To be sure, I also bounded NPOPOOA in the est file. 

Here are the relevant lines from the est file:

1  NPOPOOA     unif     10   10000      bounded output

0  RATIO_OOA_EA = NPOPOOA/2000000 hide
0  RTEA = log(RATIO_OOA_EA)    hide
0  R1 = RTEA/TDIV             hide

The 3PopExpBot20Mb.tpl file then sets growth rate to 0 at TDIV:

TDIV 1 1 0 1 0 1   //Set growth rate to 0 in deme 1
TDIV 2 2 0 1 0 1   //Set growth rate to 0 in deme 2

With these combinations of settings I do not see how population size could be going to zero (or infinity) due to negative (or positive) exponential growth. Nevertheless, the example consistently crashes  on a subset of runs with the error "free(): invalid next size (fast)".

Ivan suggested killing the demes explicitly. This should not be necessary in 3PopExpBot20Mb because there is no migration involving demes 1 and 2 earlier than TDIV. Nevertheless, to be sure, I added the following two kill lines as events to the tpl file:

TDIV 1 1 0 0 0 1   //kill deme 1
TDIV 2 2 0 0 0 1   //kill deme 2

I then ran the example 50 times from different starting parameters and a small subset of the runs crashed with the "free(): invalid next size (fast)". Is this the expected behaviour of the 3PopExpBot20Mb example? In the past I have always considered a model to be poorly coded if a subset of runs generated errors. 

Best wishes,

Jason

For reference here is the tpl file with Ivan's suggested kill lines:

#######################################################################################

//Parameters for the coalescence simulation program : fastsimcoal.exe
3 samples to simulate :
//Population effective sizes (number of genes)
NPOPAF
2000000
2000000
//Samples sizes and samples age
20
20
20
//Growth rates        : negative growth implies population expansion
0
R1
R1
//Number of migration matrices : 0 implies no migration between demes
2
//Migration matrix 0
0.0000 0.0000 0.0000
0.0000 0.0000 MIG
0.0000 MIG 0.0000
//Migration matrix 1
0 0 0
0 0 0
0 0 0
//historical event: time, source, sink, migrants, new deme size, new growth rate, migration matrix index
10 historical event
TDIV 2 0 1 1 0 1

TDIV 1 0 1 1 0 1

TDIV 1 1 0 1 0 1   //Set growth rate to 0 in deme 1
TDIV 2 2 0 1 0 1   //Set growth rate to 0 in deme 2
TDIV 1 1 0 0 0 1   //kill deme 1 which should not be necessary in this example given migration is 0
TDIV 2 2 0 0 0 1   //kill deme 2 which should not be necessary in this example given migration is 0
TBOT 0 0 0 RES1 0 1
TENDBOT 0 0 0 RES2 0 1
TENDBOT 1 1 0 1 0 1
TENDBOT 2 2 0 1 0 1
//Number of independent loci [chromosome]
1 0
//Per chromosome: Number of contiguous linkage Block: a block is a set of contiguous loci
1
//per Block:data type, number of loci, per generation recombination and mutation rates and optional parameters
FREQ 1 0 2.5e-8 OUTEXP

#################################################################################################

Here is the .est file with the bounded option added

// Priors and rules file
// *********************

[PARAMETERS]
//#isInt? #name   #dist.#min  #max
//all Ns are in number of haploid individuals
1  ANCSIZE     unif     1000  100000   bounded output
1  NBOT        unif     10   2000     bounded output
1  NPOPAF      unif     1000 100000   bounded output
1  NPOPOOA     unif     10   10000      bounded output
1  TDIV        unif     10   10000   bounded output
1  TPLUSDIV    unif     10   10000   bounded hide
0  MIG         logunif  1e-5 1e-2        bounded output

[RULES]

[COMPLEX PARAMETERS]
1  TBOT = TDIV+TPLUSDIV   output
0  RATIO_OOA_EA = NPOPOOA/2000000 hide
0  RTEA = log(RATIO_OOA_EA)    hide
0  R1 = RTEA/TDIV             hide
1  TENDBOT = TBOT+500     hide
0  RES1 = NBOT/NPOPAF     hide
0  RES2 = ANCSIZE/NBOT    hide

#################################################################################################

[J Weir]

Hi Ivan,

I am curious if you tried running these models say 50 or 100 times without the model crashing? My experience (with 3PopExpBot20Mb.tpl and my own models) is that they always crash on a at least a small subset of runs even if I actively killed them as you suggest (see my comment below where I document this). If you only ran the model once or twice then chances are you would be fine. If you did initiate many runs without them failing could you post your tpl file so I can see if I am implementing the kill step the same way as you?

Best wishes,

Jason

[Ivan Scotti]

Dear Jason,

good point. After finding where the problem was coming from, we needed to assess the variability of parameter estimation, so yes, we actually ran the simulations 100+ times. None of the runs got stuck.

Cheers,

Ivan

[Ivan Scotti]

(I'll ask my grad student Andrea to upload the last version of he TPL to avoid mistakes)

[Andrea Modica]

Dear Jason,

here's the version of the TPL we used to run a batch of 100 simulations. The script used to launch them is a costumization of the one proposed at pag 68 of the manual.

I hope it'll help you.

Cheers,

Andrea

[Laurent Excoffier]

Dear all,

beware, I have identified a bug in ver 2.7x  causing crash in simulations including exponential growth and migration. In this case, even there is no crash, it might be that simulations would be inaccurate.

An updated version will be available soon. Drop me a line if you want to test the new version before hand (win and linux version available only for the moment)

best

laurent

[Ivan Scotti]

Dear Laurent,

thank you for the warning!
Indeed our populations go through exp growth with migration (even though the combination of the two lasts only a handful of generations, and then only one pop is left).
Please let us know when we can test the new version (so we'll be able to compare with the 100 simulations we've aready run)

Regards,

Ivan

[J Weir]

Hi Laurent,

Thanks for looking into this and correcting the bug.

Cheers,

Jason

[Ivan Scotti]

Just finished running a new batch of 100 sims with fsc2705. We're now comparing the results with what we ha obtained with 2702. 

Any suggestion about where to look for potential differences in particular? Or about how to best compare the outputs - besides comparing means, shapes of parameter distributions and the like?
have a nice week

Ivan