Tskit: SLiM sometimes leaves NaNs in individuals/location

[petrelharp]

So currently SLiM does not initialize memory for individual location, and so if spatial locations are not set then these might contain NaNs. I can't provide a MWE because this depends on what was in the memory, but also note it is possible to explicitly set locations to NaN in SLiM:

initialize()
{
    initializeTreeSeq();
    initializeMutationRate(1e-2);
    initializeMutationType("m1", 0.5, "f", -0.1);
    initializeGenomicElementType("g1", m1, 1.0);
    initializeGenomicElement(g1, 0, 99);
    initializeRecombinationRate(1e-2);
}

1 { 
    sim.addSubpop("p1", 10);
    p1.individuals.x = rep(NAN, 10);
    sim.treeSeqOutput("nans.trees");
}

so that we get

>>> ts = tskit.load("nans.trees")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/peter/.local/lib/python3.7/site-packages/tskit-0.2.0a4-py3.7-linux-x86_64.egg/tskit/trees.py", line 1524, in load
    return TreeSequence.load(path)
  File "/home/peter/.local/lib/python3.7/site-packages/tskit-0.2.0a4-py3.7-linux-x86_64.egg/tskit/trees.py", line 2030, in load
    ts.load(str(path))
_tskit.LibraryError: Location values must be finite numbers

I am unclear about the extent of the problem; it may be platform-dependent. This is creating problems for the workshop that @bhaller is running in a week-ish - the tree sequence provided on that page will not load; and presumably nonspatial simulations produced during the workshop will sometimes also not load. Perhaps @bhaller can provide some other indication of how easy it is for this to happen when locations are simply not initialized.

So, it looks like we need either a SLiM release that zeros out locations or a tskit release that omits the check. This check for finiteness of locations was introduced here; at the time we said that if there was a good reason to allow NaNs in location, then we'd consider it. So, I'm opening this to decide whether we want to allow NaNs in location.

Here's arguments in either direction:

1. If we disallow NaNs then future tskit code to calculate things with locations or plot based on locations do not have to guard against them.
2. If we allow NaNs then they can be used as a "missing value", which could be useful.

A particular scenario when we might want to do (2): suppose that we have a nonspatial population that later colonizes a spatial region. If locations of the nonspatial individuals are nan, then e.g. plotting locations with matplotlib will naturally omit the nonspatial individuals. If nans are not allowed, we'd have to sort out that those values don't mean anything some other way.

Another scenario is when we're inferring a tree sequence from data, and we want to use location, but it is actually missing for some individuals.

Thinking through this more, I'm coming down on the side of allowing NaNs: in general, it is useful and important to allow recording of missingness.

And: apologies for not spotting this issue earlier, or alerting @bhaller to this proposed change.

[bhaller]

Re: how SLiM ought to behave, I would note two things. One, if we want to provide a default value for location data for non-spatial simulations then it would probably be NaN (no other value makes as much sense), so from that perspective NaN ought to be allowed. Two, I really don't want to require SLiM to initialize these values; if you have 1 million individuals in a model, with x y z locations that are each 8 bytes (double), that requires zero-initializing 22 MB of memory, if I did my math right, every generation. That would introduce a considerable amount of overhead for non-spatial simulations. I quite deliberately leave these fields uninitialized; I don't think it is unreasonable to expect the user to simply not use spatial location data if the model they are analyzing is a non-spatial model. But in any case the first argument seems dispositive to me; if we ever do make SLiM initialize these values, it will probably be to NaN.

[jeromekelleher]

Quick question: why not output 0D locations rather than nans? That way there's no initialisation required and tskit is happy.

[jeromekelleher]

Locations are a ragged array, you can mix 1d with 2d with 0d if you want.

[bhaller]

Well, @petrelharp and I decided to always write out the location values, whether the model is spatial or not, back when we did this stuff in the first place; I don't recall exactly why, but presumably simplicity/consistency factored in. Then there is also the fact that (1) the x/y/z properties are available for use even in non-spatial models, and (2) the x/y/z properties do not, in fact, actually have to be used even in spatial models – there is no requirement that the model ever set or use them. So writing out no location data in non-spatial models would be perhaps making an unwarranted assumption. I don't feel incredibly strongly about that; but neither do I see a compelling argument for changing it. And in any case, the fact remains that even if such revisions were made, NaN still seems like a reasonable/permissible value for models to use for spatial locations if they wish to. I see no reason to restrict that arbitrarily; it does nothing but remove flexibility for the user.

[jeromekelleher]

Well, IMO, if bioinformatics has taught us anything it's that flexibility is not necessarily a good thing. The point of enforcing restrictions on input data is that users should be able to analyse tree sequences with the same code regardless of what program produced them. All of this is about providing users with a consistent experience. I don't think it can be a good thing that sometimes SLiM produces tree sequences that contain uninitialised garbage for locations, and there's no obvious way of knowing whether its garbage or not, as there's not necessarily going to be NaNs in there. (This is arguably a security hole: who knows what could be in this uninitialised data? I guess modern OS's zero out pages before giving them out so it's not going to contain any passwords or anything, but still.)

Whether users should be allowed explicitly insert NaNs or not is less clear cut. Since we're not actually using the location data within tskit and even if we did it doesn't violate the model integrity (the worst that's going to happen is that we get some NaN results from calculations), I think we are being over-zealous at the moment.

So, I suggest that we drop the NaN/finite check on location data, and push out a point-release of tskit ASAP. This will only take a few hours, so everything will be sorted out in time for the workshop.

@petrelharp, do you agree? A quick point release here will be helpful anyway to get the improved performance in #349 in for benchmarks.

[petrelharp]

So, I suggest that we drop the NaN/finite check on location data, and push out a point-release of tskit ASAP. This will only take a few hours, so everything will be sorted out in time for the workshop.

I agree, this is the way to go.

[jeromekelleher]

Ok, great, I'll push out the release today. Conda packages should follow quickly.

[bhaller]

Thanks for the quick turnaround!