Efficient metadata

[bhaller]

Draft fix to make metadata access efficient in recipe 18.10; see tskit-dev/tskit#3472 (comment).

@petrelharp The function that uses the metadata gets called twice in this script, so it still accesses the metadata twice. I wasn't sure whether it would make sense to define a global, or pass the value as a parameter, etc., to get the number of metadata accesses down to one; for this script it seems like overkill, but maybe it would provide a better example in terms of pedagogy – really emphasizing that users should try to get their metadata accesses down to the bare minimum. Let me know what you think the best thing is here.

Once this PR is finalized and approved I'll modify the manual's discussion of this recipe to emphasize that minimizing metadata access is important. :->

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 76.88%. Comparing base (0657f03) to head (37054fe).

Additional details and impacted files

@@              Coverage Diff               @@
##           multitrait     #658      +/-   ##
==============================================
- Coverage       76.89%   76.88%   -0.01%     
==============================================
  Files             117      117              
  Lines           78701    78701              
  Branches        14287    14283       -4     
==============================================
- Hits            60515    60509       -6     
- Misses          18186    18192       +6     

see 10 files with indirect coverage changes

🚀 New features to boost your workflow:

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

[petrelharp]

I changed the base of this PR to point to multitrait, FYI, which seems to be what you want.

[bhaller]

I changed the base of this PR to point to multitrait, FYI, which seems to be what you want.

Huh, how odd! I could swear I set that up when I made the branch – but perhaps not. :-O Thanks!

[bhaller]

Are you otherwise happy with my fix, @petrelharp?

[petrelharp]

This is not necessary any more! I see that below we recapitate, implying that this tree sequence has all the input roots, and so the heights of all the roots should all be equal to the tick value (+/- 1) already!

[petrelharp]

See my PR to this PR, #659.

[bhaller]

OK, hmm. If you look at recipe 18.10 in the draft SLiM 6 manual (or recipe 17.10 in the SLiM 5.2 manual, which I think is pretty much the same), the goal is to see a plot showing the tree height along the chromosome before and after recapitation. Post-recapitation that's easy, but pre-recapitation I guess there was some issue that required patching the tree heights when coalescence hasn't been reached. You write "A lot of that wasn't necessary any more!" Did something change, in tskit or pyslim or somewhere, that made what we were doing unnecessary? Or was it unnecessary all along? I'd just like to understand this better, for fear of accidental breakage. I'm pretty sure that you were the primary author of this recipe, back in the day – I don't think I wrote this Python code! – so was past-you incorrect, or did the ground shift in some way? (And also: is the new code in fact better than the old code? It is actually one line longer than the old code, but I don't really understand what the new code is doing.)

I'm happy to accept your PR, I just want to double-check that all of this truly makes sense and is correct. Also, please provide a new paragraph for the manual write-up that explains what your new code is doing, and why. The old text looks like this:

Then we define a function that calculates tree heights along the chromosome. This is similar to what we did in section 18.4, but here we have to be a bit smarter because of those original ancestors preserved in the tree sequence. Every tree will have a root in one of those original ancestors, but we are interested in whether the tree has coalesced below that original ancestor (and if so, at what height), or if the tree still has multiple roots, indicating that it has not yet coalesced. We therefore check the number of roots for the tree, and if there is just a single root, we use the height of its child when the root has only a single child (indicating that coalescence occurred in the child, not in the original ancestor).

Thanks!

[petrelharp]

Yes, apologies. Let's see. The only real difference in the code is that I'm taking the height of the root(s) even if it's uncoalesced, instead of looking this up as the tick value. I don't remember why we did it that way before, but enough things have shifted it's not worth tracking down. The other difference is that the code is now checking that all roots (if there's more than one) are at the same time; this guards against, for instance, someone having simplified.

I also changed the code to output a single plot, with axis labels, to a PNG, maybe you don't want to do that?

And, I've checked out some examples and verified that (after too much time being confused by a bug) the code is doing the right thing (still); here's the output:

[bhaller]

Yes, apologies. Let's see. The only real difference in the code is that I'm taking the height of the root(s) even if it's uncoalesced, instead of looking this up as the tick value. I don't remember why we did it that way before, but enough things have shifted it's not worth tracking down. The other difference is that the code is now checking that all roots (if there's more than one) are at the same time; this guards against, for instance, someone having simplified.

OK. In the context of this specific recipe the check for all roots having the same time is unneeded, but I guess it's pedagogically useful, sure. And now I see what it's doing;set() constructs a set, which is effectively unique(). OK, I can write doc for that.

I also changed the code to output a single plot, with axis labels, to a PNG, maybe you don't want to do that?

Yeah, I don't want to do that, because the y-axis for the two plots is (intentionally and usefully) on different scales; with a single plot it is difficult to see what the pre-recapitation heights really look like. (In the manual the plot is combined, but that is using a y-axis scale that uses a cube root, which is a bit unconventional but brings out the needed detail. :->).

And, I've checked out some examples and verified that (after too much time being confused by a bug) the code is doing the right thing (still); here's the output: ...

Groovy. When I got a chance I'll merge this in, fix the plotting back the way it was (or you can fix your PR here, if you have a moment), and merge it in. Thanks!

[bhaller]

Merged, done, thanks!