Drop smck hull_offset defaults

[hossam26644]

closes #2443

[hossam26644]

@jeromekelleher
Please check this: the test test_gaps in test models fails for the smc_prime (now mapped to smc(k=1)). This tests overlapping segments in a non-discrete genome. I think smc(k=1) is not equavelent to the smc_prim in the case of a non-discrete genome. What do you think?

[jeromekelleher]

See comment above - that should make it less painful

[jeromekelleher]

I think we should make hull offset a required positional arg, like

def __init__(self, hull_offset, *, duration=None):
    ...

I don't think we put much thought into this when it was going in initially, just put something in to get things working.

(unless this breaks other things?)

[hossam26644]

it does not break other things 👍

[jeromekelleher]

Please check this: the test test_gaps in test models fails for the smc_prime (now mapped to smc(k=1)). This tests overlapping segments in a non-discrete genome. I think smc(k=1) is not equavelent to the smc_prim in the case of a non-discrete genome. What do you think?

Ooh - that's a nasty one! Can you confirm by making k=1.00001 and rerunning?

[hossam26644]

ts = msprime.simulate(
                sample_size=sample_size,
                recombination_rate=recombination_rate,
                random_seed=random_seed,
                model=msprime.SmcKApproxCoalescent(hull_offset=1.00001),
            )
            edgesets = sorted(ts.edgesets(), key=lambda e: (e.parent, e.left))
            num_found = 0
            for j in range(1, len(edgesets)):
                r = edgesets[j - 1]
                s = edgesets[j]
                if r.right != s.left and r.parent == s.parent:
                    num_found += 1
            assert num_found == 0 #fails

hull_offset=1.00001 gives similar results to 1. both fail the test

[jeromekelleher]

Hmm, good catch. We'll have to think about that.

[jeromekelleher]

Ah right, the problem is k is too big. Here, the genome is one unit long, so k=1 means we allow all possible overlaps. If we put k=0.00001 in this example we get the right answer.

Now, I'm not entirely sure what the right value of the smc_prime is in a continuous genome model. I guess it's a very small number? It's an annoying one.

[hossam26644]

In a math term, it should be when k -> 0. I am not sure if we set it to the smallest float that could be stored, it would be a sustainable universal solution. Also, we would have to set the k value differently for discrete genomes (k=1).

[jeromekelleher]

Let's try a small value like 1e-14 or something so as the default, and see if that does the right thing for the discrete genome case. Having a different mapping for discrete genome and not would be quite tricky, as we'd have to pass the discrete_genome parameter around a lot, and it would be confusing to document anyway.

Using very small values like flt_min or whatever could have unexpected consequences, so I think a reasonable like 1e-14 is less likely overall to bite us.

[hossam26644]

It does not break unit tests, and I am checking the stats verifications. Verifications does not look ok though. Check below the number of trees between msprime.SmcKApproxCoalescent(hull_offset=1e-14) and msprime.SmcPrimeApproxCoalescent()

for: the test_smck_vs_smckapprox_recomb_discrete_hotspots

recombination_map = msprime.RateMap(
           position=[0, 100, 500, 900, 1200, 1500, 2000],
           rate=[0.00001, 0, 0.0002, 0.00005, 0, 0.001],
       )
       self._run(
           samples=10,
           population_size=1000,
           recombination_rate=recombination_map,
           num_replicates=300,
           discrete_genome=True,
       )

What about we restrict the smc_prime model to discrete genome from the interface?

[jeromekelleher]

What about we restrict the smc_prime model to discrete genome from the interface?

I considered this, but it's not trivial because of the complex structure of the model argument and also it will likely break code.

Hmm, this is annoying.

OK, let's drop the remapping of the strings for now, and keep the reparameterisation of smck to require the parameter. At least this excecise has shown us that a default of 1 was definitely not the right thing to do...

It may be simpler to just do this in a different PR or update this one, it's up to you.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.22%. Comparing base (1526582) to head (902179d).
⚠️ Report is 2 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #2445   +/-   ##
=======================================
  Coverage   91.22%   91.22%           
=======================================
  Files          20       20           
  Lines       11847    11847           
  Branches     2300     2300           
=======================================
  Hits        10807    10807           
  Misses        569      569           
  Partials      471      471           

Flag	Coverage Δ
C	91.22% <100.00%> (ø)
python	98.63% <100.00%> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[jeromekelleher]

Looks good. Spotted one thing.

We can probably get rid of the hull_offset named argument in new code as well and just use SmcKApproxCoalescent(0) etc.

[jeromekelleher]

Should resinstate this test to check that it raises a TypeError (or whatever it is)

# mandatory positional argument
with pytest.raises(TypeError):
        msprime.SmcKApproxCoalescent()

[hossam26644]

check now. Also check the recommendation in the documentation: for smc_prime the documentation recommends smck(1) for discrete and smck(1e-14) for continuous genomes

[jeromekelleher]

Docstring is malformed, that's why the docs build is failing

[hossam26644]

yes, should be fixed now

[hossam26644]

@jeromekelleher Can be merged