unreferenced mutation from multiple reload

[petrelharp]

I suspect there is a simpler way to make this happen, but I haven't been able to arrive at it. If you first run this script to create in.trees:

initialize()
{
    setSeed(23);
    initializeSLiMOptions(nucleotideBased=T);
    initializeTreeSeq(timeUnit="generations");
	 initializeSex();
	 initializeMutationTypeNuc("m1", 0.5, "f", 0.0);
	 initializeMutationType("m2", 0.5, "n", 0.0, 0.1);
	 initializeGenomicElementType("g1", m1, 1.0, mmJukesCantor(8e-2));
    defineConstant("L", 100);
	 initializeChromosome(1, L, type="X", symbol="X");
    initializeGenomicElement(g1, 0, L-1);
    initializeAncestralNucleotides(randomNucleotides(L));
    initializeRecombinationRate(1e-2);
}

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

10 late() {
    sim.treeSeqOutput("in.trees");
}

And then you reload it into SLiM and write it out again immediately as out.trees by running this script:

initialize()
{
   setSeed(25);
   if (!exists("TREES_FILE")) defineGlobal("TREES_FILE", "out.trees"); 
   if (!exists("TREES_IN")) defineGlobal("TREES_IN", "in.trees"); 
	defineConstant("L", 1e2);
	initializeSLiMOptions(keepPedigrees=T, nucleotideBased=T);
   initializeTreeSeq(timeUnit="generations");
   initializeSex();
	initializeMutationTypeNuc("m1", 0.5, "f", 0.0);
   initializeMutationTypeNuc("m2", 0.5, "n", 0.0, 0.1);
	initializeGenomicElementType("g1", m1, 1.0, mmJukesCantor(4e-2));
   initializeChromosome(1, L, type="X", symbol="X");
	initializeAncestralNucleotides(paste0(rep("A", L)));
	initializeGenomicElement(g1, 0, L-1);
	initializeRecombinationRate(1e-2);
}

1 late() { 
    // output immediately! Nothing should have changed.
    md = treeSeqMetadata(TREES_IN, userData=F).getValue("SLiM").getValue("user_metadata");
    sim.readFromPopulationFile(TREES_IN);
    sim.treeSeqOutput(TREES_FILE, simplify=F, metadata=md);
}

2 early() {
    sim.simulationFinished();
}

And then you re-run the script above as follows: with the script named restart.slim, running

slim -d "TREES_IN='out.trees'" -d "TREES_FILE='out2.trees'" restart.slim

then I get

ERROR (Species::__TabulateMutationsFromTables): mutation id 809 referenced by a derived state is not in the mutation metadata table.

Error on script line 22, character 8:

    sim.readFromPopulationFile(TREES_IN);

Seems like some important mutational info read in from in the top-level metadata didn't make it into wherever it needed to be to written out.

[bhaller]

Ah, yes, hmm. SLiM now keeps a list of mutations that are "retained by the tree" - they are no longer segregating, but are referenced by ancestral derived states in the tree. When a .trees file is read in, mutations that are retained by the tree need to be put back into that list, so that they again get written out when the tree sequence is written out again. That is not presently happening; an oversight. The question, I guess, is how to recognize such mutations. Probably every mutation in the table that is not referenced by an extant individual can be assumed to be "referenced by the tree" and should be created and added to the list. I think that ought to be straightforward to do with a post-sweep through the table. Thanks.

[petrelharp]

A question here is whether to clean out unreferenced mutations from the metadata. The benefit to leaving mutations that are unreferenced by extant individuals is that we don't have to figure out which ones they are. The downside is that if someone's running a script where they repeatedly simulate-save-reload then simplification won't get rid of all the unused mutations (their metadata anyhow). How bad this is depends on how long the simulation steps are.

It doesn't seem bad to me to go through the derived states at save time, and get a list of all the SLiM IDs, and write out only those?

This raises the question that maybe you haven't thought of: what if there's a mutation in metadata that isn't referenced by any extant individuals, but has a mutation ID larger than the max among all the referenced ones? Then can we get a duplicate mutation ID in metadata? (I recall having to think about this before, but I don't remember what the conclusion was.)

[bhaller]

I believe there is existing code that checks for the max referenced mutation ID. I'll do a pass to (try to) check that that code remains effective.

Re: cleaning out, I can do this at save time, OK. I'm happier doing it at save time than at load time; since it is more common to save once and reload many times (e.g., conditional simulations) than to save many times, perhaps. Overhead either way, but the conditional-simulation case is the case where I feel more concerned about overhead.

The drawback to cleaning out only at save time is that if someone does, say, simulate-save and then simplifies away a bunch of stuff on the Python side and saves over there, then reloads into SLiM, the table will not be cleaned out, and each reload of that tree sequence will take longer than it needs to, and will occupy more memory than it needs to. Unless you intend to provide a way to clean out the table on the Python side, like with a pyslim method for doing a simplify-plus-clean? Obviously all this stuff is a pain. Personally I'm inclined to say (a) SLiM and Python never clean out the table automatically, and (b) there is perhaps a pyslim method that scans and cleans the table, but it is up to the user to choose to call that method if they care. To me that seems sufficient. Do we really want to do more than that?

[petrelharp]

Okay, so your point is that if someone does SLiM -> python -> SLiM -> ... then usage could grow. Fair enough. It's quite easy to add something to remove unneeded ones, I'll add that to the list.

I do have a strong feeling that the SLiM -> SLiM -> SLiM ... workflow might run into major problems; and if the only way to clear out unneeded things was via python then we're requiring people to drop into python in a workflow that doesn't otherwise use it. I guess my vote is to (a) put this in; (b) see if it makes saving noticably slower (I suspect not); and (c) if it does make it optional I guess? And (d) if when looking at how to clean things out, it seems like a giant pain, then revisit (but I think it'll be straightforward: just a single pass through the derived states and then through the mutations)?

[bhaller]

How is the SLiM -> SLiM -> SLiM ... workflow a problem? (By that you mean repeatedly saving from SLiM and then reloading back into SLiM, sequentially, right?) SLiM will only save out mutations into the table if they are retained by the tree; and if they're retained by the tree at save time, they should continue to be retained by the tree at load time. I'm missing where the problem comes in there. (OK, sometimes mutations that SLiM considers "retained by the tree" can in fact be discarded, because the branch they're on has been simplified away entirely, and there are not even any ancestral nodes left that reference those mutations. But if that is the problem you're concerned about, that comes up even for a pure-SLiM design that never saves/reloads at all, I think. To fix that, we'd have to do a post-pass after every simplify that scans through all derived states, and does a mark-and-sweep of the mutation IDs in use to get rid of those that are no longer referenced. We can do that, but it'd be costly, and I'm really not convinced that this is a big problem for typical models.)

AFAICT, the problem only comes up when the user goes to Python and does something that renders some of the mutations in the table unnecessary – in particular, simplifying away a chunk of the tree sequence. Am I missing the point? Maybe I got the point, above, in fact, and have now lost my grip on the point again here, lol. Sorry, my brain is mostly in issue #641 right now, lol; that is a tricky one!