Msprime: Problems with XiVsHudsonTest in verification.py

[jeromekelleher]

The refactoring of the verification tests in #1086 highlighted some issues in the tests for the Xi coalescent comparisons with Hudson here.

Basically it looks like there's a time scale problem here?

@JereKoskela, @TPPSellinger can one of you look into this please?

Please see the top of the verification.py file for instructions on how to use the new verification.py interface (please let me know if anything is unclear!)

[TPPSellinger]

This could be a remnant from when recombination and mutation where scaled differently (which I think they now are scaled the same). I will go through verification.py and let you know if this is the case.

[JereKoskela]

I'd never noticed this class of tests even existed! I think I know what's going on though.

In the previous implementation (back when there was a reference_size) we were simulating recombination in coalescent units and mutation in generations, and had to multiply one of the rates by N to account for the discrepancy. That discrepancy is gone, but it was corrected in a way that makes both trees and sites scale like N^2, rather than N. If that's right, then comparing a DiracCoalescent with population size N to a Kingman coalescent with population size N^2 should be the right thing to do.

@TPPSellinger, I'll leave this one to you if you're happy to have a look.

[TPPSellinger]

I think you are right. I will modify the verification.py, check this and let you know before tonight.

[jeromekelleher]

Great, thanks @TPPSellinger - no hurry at all here.

[TPPSellinger]

Yes I might need until tomorrow.

[jeromekelleher]

Yes I might need until tomorrow.

or next week is fine, seriously, there's no hurry!

[TPPSellinger]

Okay, I think I know what the error is here. Like @JereKoskela said, we should compare Hudson (Ne=N²) and Xi_dirac (Ne=N). However, this doesn't seem to solve the problem. I still have a factor 4 issue. Could it be that the branches of the dirac are 4 times to small (since everything seems happenning 4 times less) ?

[JereKoskela]

You're right @TPPSellinger! msp_get_common_ancestor_waiting_time_from_rate has an extra factor of 4 compared to msp_dirac_get_common_ancestor_waiting_time_from_rate. So the right comparison should be a Dirac-coalescent with population size N and a Kingman coalescent with population size N^2 / 4.

[TPPSellinger]

Okay, I fixed it. It took me some time to understand that even if the scaling error is 4, the Ne should be scaled by 2. Here is the corrected tests. Should I do a PR or it's easier for you to correct it @jeromekelleher ?
~~~~
class XiVsHudsonTest(Test):
"""
Test that Xi coalescents are equivalent to the Hudson model in the
appropriate regimes.
"""

def _run(self, xi_model, **kwargs):
    df = pd.DataFrame()
    for model in ["hudson", xi_model]:
        simulate_args = dict(kwargs)
        simulate_args["model"] = model
        model_str = "hudson"
        if model != "hudson":
            model_str = "Xi"
            simulate_args["Ne"] = 2*(math.sqrt(int(simulate_args["Ne"])))
        logging.debug(f"Running: {simulate_args}")
        replicates = msprime.simulate(**simulate_args)
        data = collections.defaultdict(list)
        for ts in replicates:
            t_mrca = np.zeros(ts.num_trees)
            for tree in ts.trees():
                t_mrca[tree.index] = tree.time(tree.root)
            data["tmrca_mean"].append(np.mean(t_mrca))
            data["num_trees"].append(ts.num_trees)
            data["num_nodes"].append(ts.num_nodes)
            data["num_edges"].append(ts.num_edges)
            data["model"].append(model_str)
        df = df.append(pd.DataFrame(data))

    df_hudson = df[df.model == "hudson"]
    df_xi = df[df.model == "Xi"]
    for stat in ["tmrca_mean", "num_trees", "num_nodes", "num_edges"]:
        v1 = df_hudson[stat]
        v2 = df_xi[stat]
        sm.graphics.qqplot(v1)
        sm.qqplot_2samples(v1, v2, line="45")
        f = self.output_dir / f"{stat}.png"
        pyplot.savefig(f, dpi=72)
        pyplot.close("all")

def test_xi_dirac_vs_hudson_recombination(self):
    self._run(
        msprime.DiracCoalescent(psi=0.99, c=0),
        sample_size=50,
        Ne=10000,
        num_replicates=1000,
        recombination_rate=0.001,
    )

def test_xi_dirac_vs_hudson_single_locus(self):
    self._run(
        msprime.DiracCoalescent(psi=0.99, c=0),
        sample_size=10,
        Ne=10000,
        num_replicates=5000,
    )

~~~~

[jeromekelleher]

Would you mind making a PR please @TPPSellinger? You might as well get git-credit for fixing the problem!

[TPPSellinger]

Sure! But it will have to wait tomorrow morning if it's okay (btw what is git credit?). May be before I do the PR I will go through all test linked to Multiple Merger and check everything is fine (or you already did that ?).

PS: I really like the new verification.py workflow

[JereKoskela]

May be before I do the PR I will go through all test linked to Multiple Merger and check everything is fine (or you already did that ?).

Thanks, that's a good idea if you have the time and energy! I suggest skipping the XiGrowth tests though. I'm still supposed to be fixing those as per #1091.

[jeromekelleher]

Thanks @TPPSellinger, you're super helpful!

Git-credit is just my shorthand for making sure that your contributions are reflected in the contributors page.

[TPPSellinger]

All tests look fine. Some of the outputs are ".png.png" but I will fix this and then do the PR. By the way, there is a "Fix Me" for the XiGrowthTest while the results looks okay (although no label).

[TPPSellinger]

Should I solve the "FIXME" by adding the x-y label and describe what it's doing before the PR @JereKoskela (or you want to do it ) ?

[jeromekelleher]

I think the XiGrowth tests have been fixed up now @TPPSellinger. If you update your branch you should see the changes that were made in #1096

[TPPSellinger]

My bad.

[TPPSellinger]

Let me know @jeromekelleher If I managed to do the PR correctly.

[TPPSellinger]

Should we close this issue @jeromekelleher ?

[JereKoskela]

I think it's safe to close this now.

Closed in #1098.