Msprime: [Sweeps] Assertion failures

[jeromekelleher]

This code fragment leads to an assertion failure:

import msprime
import numpy as np

refsize = 1e4
seqlen = 1e4
mu = 2.5e-4
rho = 2.5e-4

mod = msprime.SweepGenicSelection(
    reference_size=refsize,
    position = np.floor(seqlen/2),
    start_frequency=1.0 / (2 * refsize),
    end_frequency=1.0 - (1.0 / (2 * refsize)),
    alpha=1000.0,
    dt=1.0/(40 * refsize)
)
ts = msprime.simulate(
    10,
    model=mod,
    length=seqlen,
    recombination_rate=rho,
    mutation_rate=mu,
    random_seed=1,
    # Change to Hudson after sweep finishes
    demographic_events=[msprime.SimulationModelChange()])

We get

python3: lib/msprime.c:2127: msp_recombination_event: Assertion `x != NULL' failed.

I also got

python3: lib/msprime.c:4299: msp_change_label: Assertion `node != NULL' failed.

on a different seed. Looks like some basic bug in model changing.

[jeromekelleher]

I think there's two things happening here. The first is that we've got crazy high levels of recombination during the sweep phase, so we must be missing a scaling factor here somewhere. We're getting lots of very short segments.

We shouldn't bomb out on an assertion even in this case though, so I'd like to get to the bottom of what's happening here and return a proper error. After doing a bit of digging, it looks to me like we're hitting some limit of resolution on the msp_init_segments_and_compute_breakpoint which we should try to detect and understand.

Would you mind taking a look @daniel-goldstein? Also, if you're in there, I think the rejection sampling should be done with a simpler do-while loop. I'm not convinced the current recursive method is doing what it should.

[daniel-goldstein]

The resolution problem was definitely a top-level issue, but a symptom of these segments getting very short. The rejection sample condition tested k == y->left and we were getting k that printed the same but was slightly smaller. It is probably safer to reject if k <= y->left && y->prev == NULL (in addition to using a do-while).

I haven't yet identified the deeper problem of why we're getting such short segments. However, I am suspicious of the code in msp_run_sweep that deals with recombination mass.
https://github.com/tskit-dev/msprime/blob/c9a27e9f825921dfe586810736af01dc7813833b/lib/msprime.c#L4392-L4419

I'm also confused about the units of e_sum. We initially say that p_rec_b = rec_rates[0] * sweep_dt, leading me to think that rec_rates does not have the same units as the p_* variables. Later, we initialize e_sum = p_coal_b, but then add e_sum += rec_rates[0]:

https://github.com/tskit-dev/msprime/blob/c9a27e9f825921dfe586810736af01dc7813833b/lib/msprime.c#L4425-L4445

Can you tell if there's something fishy here in the converting of rates, or explain what's happening here? I can never remember exactly how time and recombination mass are related in terms of the rates.

[daniel-goldstein]

https://github.com/tskit-dev/msprime/blob/c9a27e9f825921dfe586810736af01dc7813833b/lib/msprime.c#L4418

also worries me, although I don't think it has to do with the problem at hand. For the couple minutes I ran the above sim, total_rate stayed between 0.000064 and 0.000117.

[jeromekelleher]

Great, thanks @daniel-goldstein. I think we should fix up the rejection sampling issue first - I think there's few scaling problems that need to be shaken out here. I'd like to see a proper error returned here.

[jeromekelleher]

Also, how did you do those lovely code quotes? They're great!

[jeromekelleher]

@andrewkern, any thoughts on the issues @daniel-goldstein points out above re scaling?

[daniel-goldstein]

Great, thanks @daniel-goldstein. I think we should fix up the rejection sampling issue first - I think there's few scaling problems that need to be shaken out here. I'd like to see a proper error returned here.

So properly rejection sample and replace the assert that was failing with say a MSP_ERR_BAD_STATE?

[daniel-goldstein]

Also, how did you do those lovely code quotes? They're great!

By highlighting line numbers (shift-click) in github files you can get a permalink and then just paste it in a comment! It's pretty great.

[jeromekelleher]

So properly rejection sample and replace the assert that was failing with say a MSP_ERR_BAD_STATE?

Yes, properly rejection sample, and return an error, but try to catch it earlier. Why is x NULL here? That wasn't clear to me.

[daniel-goldstein]

Yes, properly rejection sample, and return an error, but try to catch it earlier. Why is x NULL here? That wasn't clear to me.

x is NULL when the segment we sample, y, is at the start of a segment chain. In that case, the interval subtended by y is [y->left, y->right), which is where k has to fall. However, we can accept a smaller k when x is not NULL because then it falls in the gap still subtended by y, since y subtends (x->right, y->right). Up to precision error, k _should_ be guaranteed to fall within these intervals.

[andrewkern]

The resolution problem was definitely a top-level issue, but a symptom of these segments getting very short. The rejection sample condition tested k == y->left and we were getting k that printed the same but was slightly smaller. It is probably safer to reject if k <= y->left && y->prev == NULL (in addition to using a do-while).

I haven't yet identified the deeper problem of why we're getting such short segments. However, I am suspicious of the code in msp_run_sweep that deals with recombination mass.
https://github.com/tskit-dev/msprime/blob/c9a27e9f825921dfe586810736af01dc7813833b/lib/msprime.c#L4392-L4419

I'm also confused about the units of e_sum. We initially say that p_rec_b = rec_rates[0] * sweep_dt, leading me to think that rec_rates does not have the same units as the p_* variables. Later, we initialize e_sum = p_coal_b, but then add e_sum += rec_rates[0]:

https://github.com/tskit-dev/msprime/blob/c9a27e9f825921dfe586810736af01dc7813833b/lib/msprime.c#L4425-L4445

Can you tell if there's something fishy here in the converting of rates, or explain what's happening here? I can never remember exactly how time and recombination mass are related in terms of the rates.

finally going back through this... great catch here. this is just a bug in my implementation. everything should be scaled by sweep_dt and so the p_* variables should be used in the rejection sampling. i can fix this as a patch, or include it in my next PR

[jeromekelleher]

Great - hopefully the scaling stuff is simpler now @andrewkern, as we no longer have the idea of "model time". (I can bring you up to speed on this if you missed the changes)

[andrewkern]

yeah got some time for me today? a bunch of sweep stuff broke through these changes... trying to piece it together still

[jeromekelleher]

Sure. I've sent you an email.

[jeromekelleher]

closed in #1101