Skip to content

Bug in Dirac coalescent with single lineage in a population #2380

@jeromekelleher

Description

@jeromekelleher

In #2303 we added a test that runs a simple simulation with two populations and migration for a range of models and we hit an error with the Dirac coalescent:

jk@empire$ valgrind ./build/test_ancestry  test_simulation_replicates_dirac 
==43395== Memcheck, a memory error detector
==43395== Copyright (C) 2002-2017, and GNU GPL'd, by Julian Seward et al.
==43395== Using Valgrind-3.18.1 and LibVEX; rerun with -h for copyright info
==43395== Command: ./build/test_ancestry test_simulation_replicates_dirac
==43395== 


     CUnit - A unit testing framework for C - Version 2.1-3
     http://cunit.sourceforge.net/


Suite: msprime
  Test: test_simulation_replicates_dirac ...gsl: gamma.c:1587: ERROR: domain error
Default GSL error handler invoked.
==43395== 
==43395== Process terminating with default action of signal 6 (SIGABRT)
==43395==    at 0x4CC59FC: __pthread_kill_implementation (pthread_kill.c:44)
==43395==    by 0x4CC59FC: __pthread_kill_internal (pthread_kill.c:78)
==43395==    by 0x4CC59FC: pthread_kill@@GLIBC_2.34 (pthread_kill.c:89)
==43395==    by 0x4C71475: raise (raise.c:26)
==43395==    by 0x4C577F2: abort (abort.c:79)
==43395==    by 0x49E59E0: gsl_error (in /usr/lib/x86_64-linux-gnu/libgsl.so.27.0.0)
==43395==    by 0x4AF8E93: gsl_sf_choose_e (in /usr/lib/x86_64-linux-gnu/libgsl.so.27.0.0)
==43395==    by 0x4AF8F9F: gsl_sf_choose (in /usr/lib/x86_64-linux-gnu/libgsl.so.27.0.0)
==43395==    by 0x140A82: msp_dirac_get_common_ancestor_waiting_time (msprime.c:7822)
==43395==    by 0x139983: msp_run_coalescent (msprime.c:5208)
==43395==    by 0x13C6B6: msp_run (msprime.c:6156)
==43395==    by 0x120F4D: verify_simulation_replicates (test_ancestry.c:3011)
==43395==    by 0x121314: test_simulation_replicates_dirac (test_ancestry.c:3070)
==43395==    by 0x487C10A: ??? (in /usr/lib/x86_64-linux-gnu/libcunit.so.1.0.1)
==43395== 
==43395== HEAP SUMMARY:
==43395==     in use at exit: 125,062 bytes in 415 blocks
==43395==   total heap usage: 920 allocs, 505 frees, 390,226 bytes allocated
==43395== 
==43395== LEAK SUMMARY:
==43395==    definitely lost: 0 bytes in 0 blocks
==43395==    indirectly lost: 0 bytes in 0 blocks
==43395==      possibly lost: 0 bytes in 0 blocks
==43395==    still reachable: 125,062 bytes in 415 blocks
==43395==         suppressed: 0 bytes in 0 blocks
==43395== Rerun with --leak-check=full to see details of leaked memory
==43395== 
==43395== For lists of detected and suppressed errors, rerun with: -s
==43395== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0)
Aborted (core dumped)

This appears to trace to gsl_sf_choose being called with an argument of 1.

A quick attempt to reproduce using the Python API failed, so I'm shelving this for now. Seems like an oversight as there's no real unit testing of the Dirac or Beta coalescent with population structure.

Metadata

Metadata

Assignees

No one assigned

    Labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions