Segmentation fault about msmc

[stschiff]

Have you downloaded the most recent version? I pushed an update 10 day ago.
Stephan

On 25 Sep 2015, at 12:57, juyouhui [email protected] wrote:

Hi, stschiff
When I prepare the input file for MSMC, it occurs Segmentation fault. Need I recompile on my computer.

read 1417 SNPs from file pre_input/preinput
estimating scaled mutation rate: 0.00287284
input files: ["pre_input/preinput"]
maxIterations: 20
mutationRate: 0.00287284
recombinationRate: 0.00071821
subpopLabels: [0, 0, 0, 0]
timeSegmentPattern: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
nrThreads: 64
nrTtotSegments: 40
verbose: false
outFilePrefix: my_msmc_output
naiveImplementation: false
hmmStrideWidth: 1000
fixedPopSize: false
fixedRecombination: false
initialLambdaVec: []
directedEmissions: false
skipAmbiguous: false
indices: [0, 1, 2, 3]
logging information written to my_msmc_output.log
loop information written to my_msmc_output.loop.txt
final results written to my_msmc_output.final.txt

[1/1] estimating total branchlengthsSegmentation fault

—
Reply to this email directly or view it on GitHub #16.

[grahamgower]

I have observed this bug too. I'm on 6934da6, built using DMD64 D Compiler v2.071.0. There appears to be some kind of stack corruption. I tried some older versions of msmc and have determined that the bug disappears when adding the -release compiler flag. I also tried -boundscheck=on and -boundscheck=off: both crash with -release specified and neither crash without -release.

$ gdb build/msmc
GNU gdb (GDB) 7.11.1
Copyright (C) 2016 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-slackware-linux".
Type "show configuration" for configuration details.
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>.
Find the GDB manual and other documentation resources online at:
<http://www.gnu.org/software/gdb/documentation/>.
For help, type "help".
Type "apropos word" to search for commands related to "word"...
Reading symbols from build/msmc...done.
(gdb) run -t1 --fixedRecombination GL834414.1.txt
Starting program: /home/grg/src/msmc/build/msmc -t1 --fixedRecombination GL834414.1.txt
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib64/libthread_db.so.1".
read 21 SNPs from file GL834414.1.txt
estimating scaled mutation rate: 8.25503e-05
input files:         ["GL834414.1.txt"]
maxIterations:       20
mutationRate:        8.25503e-05
recombinationRate:   2.06376e-05
subpopLabels:        [0, 0]
timeSegmentPattern:  [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
nrThreads:           1
nrTtotSegments:      40
verbose:             false
outFilePrefix:       
naiveImplementation: false
hmmStrideWidth:      1000
fixedPopSize:        false
fixedRecombination:  true
initialLambdaVec:    []
directedEmissions:   false
skipAmbiguous:       false
indices:             [0, 1]
logging information written to .log
loop information written to .loop.txt
final results written to .final.txt
[New Thread 0x7ffff6841700 (LWP 40647)]
[1/20] Baumwelch iteration

Thread 2 "msmc" received signal SIGUSR1, User defined signal 1.
[Switching to Thread 0x7ffff6841700 (LWP 40647)]
0x00007ffff754436f in pthread_cond_wait@@GLIBC_2.3.2 () from /lib64/libpthread.so.0
(gdb) c
Continuing.

Thread 2 "msmc" received signal SIGUSR2, User defined signal 2.
0x00007ffff68777f6 in sigsuspend () from /lib64/libc.so.6
(gdb) c
Continuing.
  * [1/1] Expectation Step
Thread 1 "msmc" received signal SIGSEGV, Segmentation fault.
[Switching to Thread 0x7ffff7fb8c40 (LWP 40643)]
0x00000000005a0dc0 in model.propagation_core.PropagationCore.propagateMultiForward(const(model.stateVec.State_t), model.stateVec.State_t, const(model.data.SegSite_t), const(model.data.SegSite_t)) const.__require() (this=0x7fffffffbd10)
    at model/propagation_core.d:44
44          assert(to_segsite.pos > from_segsite.pos);
(gdb) bt
#0  0x00000000005a0dc0 in model.propagation_core.PropagationCore.propagateMultiForward(const(model.stateVec.State_t), model.stateVec.State_t, const(model.data.SegSite_t), const(model.data.SegSite_t)) const.__require() (
    this=0x7fffffffbd10) at model/propagation_core.d:44
#1  0x000000000058398b in model.propagation_core_fastImpl.PropagationCoreFast.propagateMultiForward(const(model.stateVec.State_t), model.stateVec.State_t, const(model.data.SegSite_t), const(model.data.SegSite_t)) const (
    this=0x7ffff7eb8700, to_segsite=0x7ffff7f0aec0, from_segsite=0x7ffff7ebc000, to=0x7ffff7f0ae80, from=0x7ffff7f023c0) at model/propagation_core_fastImpl.d:270
#2  0x000000000059ee0a in model.msmc_hmm.MSMC_hmm.runForward() (this=0x7ffff7eb8100) at model/msmc_hmm.d:138
#3  0x00000000005a8ada in expectation_step.singleChromosomeExpectation(const(model.data.SegSite_t[]), ulong, const(model.propagation_core.PropagationCore)) (__HID138=0x7fffffffbff0, propagationCore=0x7fffffffbff0, hmmStrideWidth=1000, 
    data=...) at expectation_step.d:77
#4  0x00000000005a894e in expectation_step.getExpectation(const(model.data.SegSite_t[][]), model.msmc_model.MSMCmodel, ulong, ulong, bool).__foreachbody6(ref const(model.data.SegSite_t[])) (this=0x7fffffffd5c0, __applyArg0=...)
    at expectation_step.d:59
#5  0x000000000057fc79 in std.parallelism.ParallelForeach!(const(model.data.SegSite_t[])[]).ParallelForeach.opApply(scope int(ref const(model.data.SegSite_t[])) delegate).doIt() (this=0x7fffffffd480)
    at /usr/include/dmd/phobos/std/parallelism.d-mixin-3785:3829
#6  0x00000000005ca0d8 in std.parallelism.run!(void() delegate).run(void() delegate) ()
#7  0x00000000005c9c18 in std.parallelism.Task!(std.parallelism.run, void() delegate).Task.impl(void*) ()
#8  0x00000000005fa2b3 in std.parallelism.AbstractTask.job() ()
#9  0x00000000005c98db in std.parallelism.submitAndExecute(std.parallelism.TaskPool, scope void() delegate) ()
#10 0x000000000057fb85 in std.parallelism.ParallelForeach!(const(model.data.SegSite_t[])[]).ParallelForeach.opApply(scope int(ref const(model.data.SegSite_t[])) delegate) (this=..., dg=...)
    at /usr/include/dmd/phobos/std/parallelism.d-mixin-3785:3835
#11 0x00000000005a8844 in expectation_step.getExpectation(const(model.data.SegSite_t[][]), model.msmc_model.MSMCmodel, ulong, ulong, bool) (__HID137=0x7fffffffd730, naiveImplementation=false, maxDistance=1000, hmmStrideWidth=1000, 
    msmc=0x7ffff7ec3940, inputData=...) at expectation_step.d:57
#12 0x00000000005aa737 in msmc.run() () at msmc.d:263
#13 0x00000000005a8c35 in D main (args=...) at msmc.d:89
(gdb) p to
Internal error: `this' is not an aggregate
(gdb) up
#1  0x000000000058398b in model.propagation_core_fastImpl.PropagationCoreFast.propagateMultiForward(const(model.stateVec.State_t), model.stateVec.State_t, const(model.data.SegSite_t), const(model.data.SegSite_t)) const (
    this=0x7ffff7eb8700, to_segsite=0x7ffff7f0aec0, from_segsite=0x7ffff7ebc000, to=0x7ffff7f0ae80, from=0x7ffff7f023c0) at model/propagation_core_fastImpl.d:270
270       override void propagateMultiForward(in State_t from, State_t to,
(gdb) p to
$1 = (model.stateVec.State_t *) 0x7ffff7f0ae80

[grahamgower]

I've tried running this through valgrind, to check for memory errors, then piping the output through ddemangle. Some memory errors can be ignored, as they are from the garbage collector (http://forum.dlang.org/thread/CAE8u=e6WMin8+xt2+8M5kHt8=[e-mail]). Its clear that there are some errors in the non-release version (vg.ddemangle.txt) that do not appear in the -release version (release.vg.ddemangle.txt). Note that I did not wait for the -release run to complete, but killed it with CTRL-C after a few mins when it had completed the first Expectation Step.

release.vg.ddemangle.txt
vg.ddemangle.txt

[grahamgower]

The attached input file results in a crash within 10 seconds on my system.
GL834414.1.txt

[stschiff]

I have now built in an option to artificially constrain the coalescence rates, see release 1.1.0, with options --loBoundLambda and --hiBoundLambda. Lambda is the scaled coalescence rate, with typical values around 1000 or so (it's on average around 2/theta, where theta is the heterozygosity), so you could try setting a lower bound of 10, and an upper bound of 1e6 or so. Perhaps that resolves the issue. If not, please let me know.