Allow `keep_unary` to be passed to simplification funcs

[hyanwong]

I can imagine that it would be very handy to be able to output in tree sequence all the intermediate nodes produced by SLiM - i.e. not just the coalescence points, but all the "unary nodes" too. This would be somewhat intermediate between "saving everybody" and the normal output. Individuals would still not be saved, of course, unless they happen to be "retained" - this is just about saving the node info (genomes) themselves.

The way I imagine this working would be to add an optional keep_unary arg to treeSeqSimplify(), which could be T, F, or NULL (default). There would also be a keep_unary option to initializeTreeSeq() that would default to F and would set the default keep_unary flag during automatic simplification. This value would be the default used when calling treeSeqSimplify() with no params.

What do you think @petrelharp ? I guess you're the one doing "remember everyone" simulations.

[bhaller]

As I mentioned in the review on #25, let's try to envision the ideal final design for this stuff. Is it a keepUnary=T option on initializeTreeSeq() and a permanent=F option on treeSeqRememberIndividuals()? In other words, are these two features completely orthogonal? Or is there some interaction between them, and some better API design? Because keepUnary=T strikes me as a kind of "remembering" or "retaining" too, but we're proposing handling it in a very different way. It could be handled instead as another type of flag passed to treeSeqRememberIndividuals(), right? A flag saying "mark these individuals as being kept even if their nodes are unary"? I don't see a clear problem with what has been proposed, but I want to make sure we're going to arrive at a final design that we won't regret later.

[hyanwong]

I see them as completely orthogonal, which is why I opened up a new issue. I guess you could think of this as a sort of "remembering", but of nodes not individuals. The distinction is very clear in tree-sequence world, and I think that the "keep_unary" function should automatically play nicely (and logically) with the permanent=F option in #25 . But let's see what @petrelharp thinks too?

[petrelharp]

Not even so much orthogonal, but copmplementary: if you do keep_unary=T and you Retain all the individuals, then you'll get the individuals corresponding to retained unary nodes as well. (However, I'd imagine that would be nearly the same as just Remembering everyone.)

[hyanwong]

Not even so much orthogonal, but complementary

Yes, bad phrasing on my part. That's what I meant by "play nicely together".

(However, I'd imagine that would be nearly the same as just Remembering everyone.)

I think it would be for short simulations. For long running ones, I suspect that there will be lots of ancestors in the distant past that are genealogical ancestors but not genetic ancestors anywhere (unless e.g. there is loads of gene conversion)

[hyanwong]

Oh, one extra wrinkle I just thought of. There's not much point to allowing treeSeqSimplify(keep_unary=T) if by default the next time that simplification is done automatically, it removes all the unary nodes anyway, because we started the simulation with initializeTreeSeq(keep_unary=F). But unless anyone has a better suggestion, I guess there's no harm in allowing it, as long as it's documented that subsequent auto-simplification could remove the unary nodes.

[bhaller]

Oh, one extra wrinkle I just thought of. There's not much point to allowing treeSeqSimplify(keep_unary=T) if by default the next time that simplification is done automatically, it removes all the unary nodes anyway, because we started the simulation with initializeTreeSeq(keep_unary=F). But unless anyone has a better suggestion, I guess there's no harm in allowing it, as long as it's documented that subsequent auto-simplification could remove the unary nodes.

Why even have it as a parameter to treeSeqSimplify()? Why not just use the flag that was given to initializeTreeSeq() for all simplifications?

[bhaller]

(However, I'd imagine that would be nearly the same as just Remembering everyone.)

I think it would be for short simulations. For long running ones, I suspect that there will be lots of ancestors in the distant past that are genealogical ancestors but not genetic ancestors anywhere (unless e.g. there is loads of gene conversion)

OK, probably I'm just confused, but. Remembering all individuals would mean that whole evolutionary branches that eventually went extinct would still be preserved in the tree sequence, right? Every individual would be remembered. Whereas retaining every individual and setting keepUnary=T would keep every individual that is genetically ancestral to the extant generation, but it would allow the evolutionary branches that went extinct to be pruned away by simplify, right? So isn't there a huge difference between the two, given that most evolutionary branches eventually go extinct? Or what am I misunderstanding?

[petrelharp]

If we keepUnary, then we'll prune out only individuals who leave behind no genome to the present day. If you're doing a long chromosome, that's nearly the same as individuals who have no pedigree ancestors in the present day. The number of pedigree descndants of an individual in a WF population acts, for a while anyhow, like a branching process with a Poisson(2) offspring distribution, whose probability of extinction is about 20%. So, I'd expect if you're doing a simulation with a long genome (like, 1Morgan long), that keepUnary would end up with about 80% as many individuals as compared to if you Remembered everyone. If you're doing a short genome though, then I agree, the number would be a lot smaller.

[bhaller]

If we keepUnary, then we'll prune out only individuals who leave behind no genome to the present day. If you're doing a long chromosome, that's nearly the same as individuals who have no pedigree ancestors in the present day.

Yes – but it's much, much smaller than the total number of individuals that have ever lived, right? Which is what just remembering every individual in every generation would do, right? I mean, this is the main point of simplification, isn't it? That over time, most evolutionary branches go completely extinct: most individuals leave no pedigree descendants over the long run. But remembering every individual in every generation keeps all of those extinct branches, and effectively turns off simplification, doesn't it? Or am I misremembering how "remembering" works?

[petrelharp]

Yes – but it's much, much smaller than the total number of individuals that have ever lived, right?

Well, this is an empirical question, but my calculation above suggests that the number of individuals with no pedigree descendants is about 80% of the total number of individuals ever.

I think you're getting mixed up by the fact that the tree sequence only has to keep track of nodes with coalescences in extant material, not all nodes that are ancestors. So the condition for appearing in a tree sequence is "I'm the MRCA of at least two samples at some locus", which is a lot stricter than "I'm an ancestor of at least one sample at some locus".

[petrelharp]

Just to be clear, I'm fine with this proposed change, but I suspect that in practice I'm going to keep on remembering everyone.

[hyanwong]

Well, this is an empirical question, but my calculation above suggests that the number of individuals with no pedigree descendants is about 80% of the total number of individuals ever.

I think it's a relatively standard result that the pedigree ancestors encompass about 80% of an entire panmictic population (I think Richard Dawkins and I were the first to write about this in a popular science book, based on Chang's calculations, around the year 2000!). But that's not the same as the percentage of genetic (not pedigree) ancestors, which are the ones we would keep if we did keep_unary. Are pedigree and genetic ancestors really the same order of magnitude as @petrelharp says? I can see that's the case as the chromosomes get very long, but I'm not sure it is in e.g. humans. Graham Coop has a nice blog post about this, ISTR, and I asked him about it (?on twitter), assuming that far enough into the past, the genetic lineages get cut up into such fine pieces that they basically pass through all the pedigree ancestors. But he said that doesn't happen - something to do with the asymptotic effect of recombination.

I guess if we implement this feature, it'll be easy to find out!