saving without simplification results in mutations without metadata

[petrelharp]

The nice logic around retaining mutations or not means that we're throwing out uneeded bits in two ways: retaining mutations, and at simplification. However, retaining is assuming that we will simplify in the future, which is optional. So this recipe:

initialize()
{
    setSeed(23);
    initializeSLiMOptions(keepPedigrees=T);
    initializeTreeSeq(simplificationRatio=INF, timeUnit="generations");
    initializeMutationRate(1e-2);
    initializeMutationType("m1", 0.5, "f", -0.1);
    initializeGenomicElementType("g1", m1, 1.0);
    initializeGenomicElement(g1, 0, 99);
    initializeRecombinationRate(1e-2);
}

1 early() { 
    sim.addSubpop("p1", 10);
}

10 late() {
    sim.treeSeqOutput("out.trees", simplify=F);
}

produces a tree sequence for which not all mutations have information, as seen as follows:

import tskit
ts = tskit.load("out.trees")

md = { x['mutation_id'] : x for x in ts.metadata['SLiM_mutation_list'] }
ds = [ int(j) for x in ts.mutations_derived_state for j in x.split(",") ]
for k in ds:
    assert k in md, f"missing {k}"

Another way to get to this error would be to (a) write out the tree sequence without simplification; (b) use pyslim to mark more individuals as samples; (c) try to read it back into SLiM, who would now have the "unknown mutation" error.

[petrelharp]

A solution to this would be to actually remove most of the "retain mutations" code, and leave it up to the code at simplify time to do the cleaning. This would be equivalent (storage-wise, at least?) of what we were doing before (writing everything down in the tree sequence tables and cleaning it later). I don't know if there's performance implications otherwise.

[bhaller]

Hmm. Well, yes. Good point.

This raises the question: has it proved to be useful, to allow people to write from SLiM without simplification? Is that a thing we actually want to continue to support? You've got some experience now with hybrid simulation workflows of various kinds. Have any of the various workflows you've seen involve saving without simplification? Or is it ever useful to you or others for other reasons, like for debugging? If we could just ditch this feature, it might make our lives simpler.

Assuming we don't want to ditch the feature, then yes, it appears to force SLiM to keep a retain on all mutation objects until such time as they are simplified away. That would be a simpler design overall, but would raise the memory usage somewhat. (As you note, it is similar to what we used to do when the mutation metadata was in the mutation table; but now we would have a retain on each Mutation object until simplify time, rather than just keeping the metadata for each mutation until simplify time. The Mutation objects are somewhat more heavyweight; each one is a separate allocated block, and each one has Dictionary as a superclass so it has Dictionary-related cruft, and so forth. Overall I think this is not a big deal; typically mutations are not the main culprit for memory footprint anyway. But it will be a higher high-water mark than it was before.)

I think this design shift would be pretty straightforward. There would still need to be some concept of being "retained by the tree", since Mutation objects normally get freed when they are lost, and we need to prevent that. So instead of slapping an extra retain on only certain mutations when they are first added to the tree, we'd need to slap a retain on all mutations when they are first added to the tree, I guess. After simplification, we'd do the pass through all derived states, and do one release (the opposite of retain) on all mutations not referenced by a derived state. We'd probably always do that sweep, since we won't have a list of mutations retained by the tree any more (and thus won't have a count of them, so won't know whether we have > 10,000 or not). This will add a little new overhead, then: retaining every mutation, and sweeping on every simplify (and probably freeing most of the mutations we previously retained). Seems like we just have to accept that overhead; I don't see a better design at the present moment.

[bhaller]

(I'm realizing, as I ponder this: perhaps you were envisioning that we would write all the mutation metadata into a buffer, one row per mutation, as we simulate along, and then clean unreferenced rows out of that buffer at simplify time. I definitely don't think we want to do that; it adds lots of headaches and overhead. I want to keep the Mutation objects themselves as the holders of the metadata, as we are doing now, and copy metadata from them into a buffer when we save, as we are doing now. So the sweep at simplify is a sweep that releases, and probably frees, Mutation objects that are no longer referenced by the tree. I think that's a much nicer design, since then we don't need to worry about when we recorded the metadata, updating the metadata if the Mutation object changes later, etc.)

[petrelharp]

Good question. I'll ponder, but in the meantime, does #646 interact with this?

[bhaller]

Maybe? :->

[bhaller]

Well, the design shift was not as straightforward as I had hoped. Most of the complication in the existing design needed to be kept, and some new complications entered in as well. Overall it is not that many lines of code, but the logic is a bit subtle, and it feels brittle.

A new wrinkle that just occurred to me is that with a multi-chrom model, with one tree sequence per chromosome, the top-level mutation metadata table will need to be different for each chromosome. So each chromosome kinda needs to keep its own separate vector of "mutations retained by the tree sequence"; or when writing out each tree sequence, the code needs to filter down the overall mutation list into just the subset of mutations that are associated with the chromosome being written out. I think it has to be the former, really. Probably this will force me to divide the mutation registry kept by the Species into separate mutation registries kept by each chromosome; and similarly, the vector of substitutions kept by the Species will have to be split into separate vectors per chromosome. I've been wondering whether those things would have to happen someday; now they have happened. :-O So this is going to be a large wave of changes, but I don't see an alternative.

Related to that: as far as I can tell, there is a big existing bug in SLiM 5, and I don't understand how it has not come up before now. The mutations table is not, of course, shared between the multiple tree sequences; each chromosome has its own mutations. BUT there is only a single list of substitution objects, kept in population_.treeseq_substitutions_map_, for the whole species. In Species::RecordNewDerivedState() when making a new derived state at a given position, we just look up all substitutions at that position and add them all to the derived state we're making – regardless of which chromosome the substitutions are associated with! So if a substitution occurs at position 17 in chromosome 5, then from then on derived states generated for position 17 in ALL chromosomes will contain a reference to that substitution in them! If there is a way that the present design avoids that, I'm not seeing it. It must've just not been hit because, for most models, it is unlikely that the same position will happen to get used twice, in two different chromosomes. :-O But wow. So I'm going to need to rearchitect all of this kind of stuff to be per-chromosome. Aargh. SLiM is getting complex enough that it is hard to hold the whole design in my head, and hard to see design errors like this that ought to have been obvious. Or maybe it is all fine, and I am just confusing myself and chasing my own tail; it's hard to tell. Aargh.

[bhaller]

I've opened #649 on those multi-chrom wrinkles, since that will be a big project in itself. For now I will ignore those wrinkles, and will get this issue fixed for single-chromosome models; then I'll return to all this in #649 to fix it up for multi-chrom.