fastsimcoal returning only 0 likelihood, not exploring space

Shannon Kieran

unread,

Sep 29, 2020, 3:00:52 PM9/29/20

to fastsimcoal

Hi Laurent and Fastsimcoal folks,

I'm having a weird issue with fastsimcoal26. I'm getting some segmentation faults in my error log, but the program continues to run. It only explores a few values, shuffling them between parameters but never varying from them regardless of starting seed. It returns a best likelihood of 0. I'm allocating 60G memory total, but 4G memory max per core. These settings seem to have worked for other folks on my cluster.

I'm using the 2D-SFS from two populations (with 31 and 35 individuals) to model an isolation-with-migration scenario. I built the 2D-SFS in angsd using 3499 shared RAD loci. I added a header for fastsimcoal using an R script from Ismail Saglam originally written by Laurent (as far as I can tell).

Here is my template file (Ne is a rough estimate to "anchor" the model, I get the same results when I replace with NPOP0/NPOP1 here and in .est file):

//Parameters for the coalescence simulation program : simcoal.exe

2 samples to simulate :

//Population effective sizes (number of genes)

75025000

190275000

//Samples sizes and samples age

31

35

//Growth rates: negative growth implies population expansion

0

//Number of migration matrices : 0 implies no migration between demes

2

//Migration matrix 0

0 MIG10

MIG01 0

//Migration matrix 1

0 0

//historical event: time, source, sink, migrants, new deme size, growth rate, migr mat index

1 historical event

TDIV 1 0 1 RESIZE 0 1

//Number of independent loci [chromosome]

3499 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 gen recomb and mut rates

FREQ 1 0 2.0e-8

And here's the .est file:

// Priors and rules file

// *********************

[PARAMETERS]

//#isInt? #name #dist.#min #max

//all N are in number of haploid individuals

1 ANCSIZE unif 100 100000 output

0 N0M10 logunif 1e-2 500 hide

0 N1M01 logunif 1e-2 500 hide

1 TDIV unif 100 50000 output

[RULES]

[COMPLEX PARAMETERS]

0 RESIZE = ANCSIZE/NPOP1 hide

0 MIG10 = N0M10/75025000 output

0 MIG01 = N1M01/190275000 output

Here is the code I use to run it:

fsc26 -t ${F1}_${x}.tpl -n100000 -m -e ${F1}_${x}.est -M -L40 -c8 -q

Here is the .simparam file output from a run:

Simulation parameters for current simulations

=============================================

Deme sizes

Deme 0 75025000

Deme 1 190275000

Sample sizes

Deme 0 31

Deme 1 35

Sample ages

Deme 0 0

Deme 1 0

Growth rates

Deme 0 0

Deme 1 0

Migration matrix

Note that migration rates have been pre-multiplied by factor TDemeCollection::CTTimeUnit (1e+06)

0.0000000 58.1413473

12.6341659 0.0000000

Migration matrix

Note that migration rates have been pre-multiplied by factor TDemeCollection::CTTimeUnit (1000000.0000000)

0.0000000 0.0000000

Historical events

Event 0

#Time : 41089

#Source : 0

#Sink : 1

#Migrants : 1.0000000

#New size : 0.6317499

#New growth rate : 0.0000000

#New migr. matrix : 1

Number of independent loci to simulate : 3499

With the same chromosomal structure

Number of linkage blocks to simulate in structure 1: 1

1 partially linked FREQ data:

recombination and mutation rates : 0.0000000000 0.0000000200

Do not output expected MAF spectrum

Do not output expected DAF spectrum

^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^

And here is a part of the brent_lhoods file to show the exploration of the space, also observable in the outfile (attached):

Param# ANCSIZE NPOP0 NPOP1 TDIV MIG10 MIG01 MaxEstLhood MaxObsLhood MAFLHood_1

0 101 35494 910 38288 9.78114e-05 2.10577e-04 0.00000

1 101 98756 910 38288 3.51545e-05 2.10577e-04 0.00000

2 101 98756 98869 38288 3.51545e-05 1.93818e-06 0.00000

3 101 98756 98869 38288 1.09360e-07 1.93818e-06 0.00004 101 98382 101 24311 0.0046550 4.5920425 0.00000

5 101 98382 101 49115 0.0046550 4.5920425 0.00000

0 98852 98382 101 49115 0.0046550 4.5920425 0.00000

1 98852 101 101 49115 4.5343248 4.5920425 0.00000

2 98852 101 98852 49115 4.5343248 0.0046918 0.00000

3 98852 101 98852 49115 1.08037e-04 0.0046918 0.00000

4 98852 101 98852 49115 1.08037e-04 1.10396e-07 0.00000

5 98852 101 98852 102 1.08037e-04 1.10396e-07 0.00000

0 101 101 98852 102 1.08037e-04 1.10396e-07 0.00000

1 101 98852 98852 102 1.10385e-07 1.10396e-07 0.00000

2 101 98852 101 102 1.10385e-07 1.08049e-04 0.00000

3 101 98852 101 102 0.0046670 1.08049e-04 0.00000

4 101 98852 101 102 0.0046670 4.5677929 0.00000

5 101 98852 101 49427 0.0046670 4.5677929 0.00000

0 98852 98852 101 49427 0.0046670 4.5677929 0.00000

1 98852 101 101 49427 4.5677867 4.5677929 0.00000

2 98852 101 98852 49427 4.5677867 0.0046670 0.00000

3 98852 101 98852 49427 1.08044e-04 0.0046670 0.00000

4 98852 101 98852 49427 1.08044e-04 1.10392e-07 0.00000

5 98852 101 98852 102 1.08044e-04 1.10392e-07 0.00000

0 101 101 98852 102 1.08044e-04 1.10392e-07 0.00000

1 101 98852 98852 102 1.10392e-07 1.10392e-07 0.00000

2 101 98852 101 102 1.10392e-07 1.08044e-04 0.00000

3 101 98852 101 102 0.0046670 1.08044e-04 0.00000

4 101 98852 101 102 0.0046670 4.5677903 0.00000

5 101 98852 101 49427 0.0046670 4.5677903 0.00000

0 98852 98852 101 49427 0.0046670 4.5677903 0.00000

1 98852 101 101 49427 4.5677903 4.5677903 0.00000

2 98852 101 98852 49427 4.5677903 0.0046670 0.00000

3 98852 101 98852 49427 1.08044e-04 0.0046670 0.00000

4 98852 101 98852 49427 1.08044e-04 1.10392e-07 0.00000

5 98852 101 98852 102 1.08044e-04 1.10392e-07 0.00000

0 101 101 98852 102 1.08044e-04 1.10392e-07 0.00000

The program generally runs for about 18 hours. The segmentation faults look like this:

srun: error: bigmem4: task 1: Segmentation fault

srun: error: bigmem4: task 2: Segmentation fault

srun: error: bigmem4: task 6: Segmentation fault

They occur throughout the run. The end of the run does say "CANCELLED,FAILED" as the status but the segfaults start pretty quickly after the run.

Is this just a memory allocation issue? Is it something with my infiles? Am I asking it to do something weird? Any guidance would be appreciated. I'm trying to figure out how to attach the .obs and .out files, if I can't I'll copy them into a reply.

Laurent Excoffier

unread,

Sep 29, 2020, 3:17:35 PM9/29/20

to fastsimcoal

Hi

it seems that the instruction

//Number of independent loci [chromosome]
3499 0

is wrong. It should be

//Number of independent loci [chromosome]

1 0

Try to see if this works better

laurent

Shannon Kieran

unread,

Oct 5, 2020, 3:20:47 PM10/5/20

to fastsimcoal

Hi Laurent,

Unfortunately this did not fix the issue. The model explored the same set of numbers and I still got a segmentation fault.

In most ways this run was identical to previous runs, but the simparam file in this case suggested the size went to infinite:

Historical events

Event 0

#Time : 20491

#Source : 1

#Sink : 0

#Migrants : 1.0000000

#New size : inf

#New growth rate : 0.0000000

#New migr. matrix : 1

I can send the .obs and outfile over email. (Sorry, thought this message posted days ago!)

Reply all

Reply to author

Forward