Msprime: worry about floating point and probabilities in mutation generation

[petrelharp]

For instance, in mutation_matrix_choose_root_state:

    while (u > params.root_distribution[j]) {
        u -= params.root_distribution[j];
        j++;
        assert(j < params.num_alleles);
    }

that assert might be hit if the root distribution doesn't sum to exactly 1. Since we do error checking that the probabilities make sense, why don't we instead change this to

    while (u > params.root_distribution[j] && j + 1 < params.num_alleles) {
        u -= params.root_distribution[j];
        j++;
    }

ie, taking the last entry? (if I did that right)

[jeromekelleher]

All we want to do is make sure that we don't fall off the edge of the array, so we could do

while (j < params.num_alleles && u >  params.root_distribution[j]) {
    /* etc */
}

Short circuit behaviour guarantees that the second part of the && doesn't get evaluated.

We should have this kind of paranoia in there though.

[petrelharp]

I don't think that will do it, since the loop could still end with j equalling num_alleles, and the next line is

    site->ancestral_state = params.alleles[j];

[jeromekelleher]

Good point. We should probably just use the GSL functions for doing this anyway, so that we don't have to worry about these details.

[jeromekelleher]

Converting this code to use the GSL methods would be a good place for someone to start out on the C code - any takers?

[castedo]

I'm ready to jump on this, if nobody else wants to take. @jeromekelleher go ahead and assign it to me and I'll do it next week. @nspope were you hoping to do this? I don't want to deprive anybody of the fun! ;-)

[nspope]

Go ahead!

[jeromekelleher]

Great, thanks @castedo!

[castedo]

Good point. We should probably just use the GSL functions for doing this anyway, so that we don't have to worry about these details.

@jeromekelleher Here's an update on a slight change of course I'm currently planning. I'm planning to not use that GSL function. I'm planning to move some of this suspect code into separate functions with a calling interface similar to the GSL ones you mention. And make sure the rounding issues are handled. Maybe I'll write a unit test too.

Using that GSL function has a space-time trade-off. A new allocation gets made so that discrete random picks can be faster. This new allocation will make the code and data structures more complicated. I'm profiling test_mutations test_single_tree_mutgen_keep_sites_many_mutations with kcachegrind and I see not much time spent in the code picking random discrete choices using the matrix probabilities. Let me know if there is a better example run to profile. Or if you have a good sense that in real-world scenarios this random picking does get hit a significant amount (relative to all the other code running).

[jeromekelleher]

We're not so worried about about performance here @castedo so much as correctness. This is the main motivation for using GSL's functions - they're battle hardened, and we don't have to worry about what happens in tricky corner cases. So, unless it makes things significantly slower, I would prefer to use GSL's methods.

[castedo]

I'd be shocked if the extra GSL struct made things signif slower. I think I see good way to keep the code simple and include the GSL struct. It might make the print function output uglier, but it will be a sacrifice for more cleaner correct-code reuse.

[castedo]

In the process of debugging subtle numerical issues I ended up creating two
solutions in my heavily #ifdef'd https://github.com/castedo/msprime/tree/ifdef_party branch.

One solution is essentially equivalent to the code that @petrelharp originally proposed. This solution is cleaned up in PR #1156 and does not use that GSL function with lookup table. I have re-written Peter's orig proposed code to be a little bit clearer and easier to understand (I think). It is the simple_ran_discrete function in mutgen.c of PR #1156
https://github.com/tskit-dev/msprime/pull/1156/files#diff-2855ce2e7aad0c9d8740b3c94e165a5aR200

The two places in mutgen.c where there could be a rare rounding issue now call this function which does not have that rounding issue. I can go into more detail on how this function behaves in all the various corner cases of doubles that get passed in. Overall, I would summarize the function as perfectly good and properly behaving for this application. Or in other words, correct. This PR #1156 is passing all tests. except for a code coverage drop.

On the other hand, there is PR #1155 which does use the GSL function and lookup table. However there are a number of issues. The biggest issue is that unit tests are failing, specifically test_flip_flop_mutation and test_uniform_mutations. These tests are failing because they assume the way the random seed and number generation is used is like in the simple way the currently implementation does. The simple implmentation is re-coded in Python inside tests/test_mutation.py
PythonMutationMatrixModel.choose_allele.

The other issue is that the GSL lookup table used by gsl_ran_discrete does not preserve the probabilities passed in close to double epsilon precision (around 1e-16). For the test cases of test_BLOSUM62 and test_PAM some probabilities actually get distorted by around 1e-14. This is what caused some confusing results yesterday.

Depending on what you consider correctness, I would describe the advantage of the GSL gsl_ran_discrete function and look up table as a performance benefit, not a correctness benefit. The implementation I imagine is rather complicated. I'm not confident it would be easy to understand why it returns slightly different probabilities than it is given. I'm sure it's a great tool for getting a performance boost. The distortion probably stays very small and is good enough for the vast majority of simulations.

In contrast, simple_ran_discrete I don't see as having any downside in correctness. It's so simple it's practically a formal definition of _A_ correct behavior (there isn't just one correct behavior).

It seems to me using the GSL functions is not worth it. It doesn't really gain any benefit other than making it easy to assume a black box works correctly. In contrast simple_ran_discrete looks perfectly correct to me and is a simpler solution. So I'd go with the simple solution (and it's less work too).

[petrelharp]

This makes sense to me, @castedo! And, you're right - the GSL implementation is complex because it's designed to work well at the large-N case, which we are nowhere near. The only corner cases besides the doesn't-quite-sum-to-one one, which you've dealt with nicely, would have to do with very small probabilities (eg what if probs are 0.5 - 0.5e-12, 1e-12, 0.5-0.5e-12). But (a) there is no such model in popgen, and (b) we don't even know if GSL deals with that better than our implementation.

1156 looks good - I'll have a few comments if we decide to go that direction.

[castedo]

To be more easily understood and predictable for very small probabilities, I can suggest coding a running sum of the probabilities rather than a running subtraction from the random number. I agree these very small probabilities aren't much of a practical concern in the popgen application. But here it is:

static size_t
simple_ran_discrete(gsl_rng *rng, size_t num_probs, double const *probs)
{
    double u = gsl_ran_flat(rng, 0.0, 1.0);
    size_t i = 0;
    double sum = 0.0;
    while (i + 1 < num_probs) {  /* ignore last probability, don't care */
        assert(probs[i] >= 0.0);
        sum += probs[i];
        if (sum > u)
            break;
        ++i;
    }
    return i;
}

The reason is that the probabilities are fixed and predictable. This makes the running sum and associated IEEE roundings easier to imagine as there are only a few predictable calculations. In contrast, by keeping a running subtraction of the random number, there is a random and crazy large number of possible IEEE roundings that can occur. And this makes it much more difficult to think about the cases with very small probabilities.

In this implementation it's easy to see that if a probability is very small [1] then its index will never be returned. To make this truly hold I also changed the non-strict comparison to a strict comparison for the break ... to cover the crazy rare corner case of gsl_ran_flat returning zero when the first probability is also zero :).

[1] "very small" in this case meaning the number is so small that it causes no change in the cumulative sum of the user provided probabilities. Numbers greater than DBL_EPSILON will not be this "very small" and numbers less than DBL_EPSILON might be depending on the running sum.

The chances I'm wrong about something here is greater than gsl_ran_flat returning zero and DBL_EPSILON. ;-)

[petrelharp]

Hey, I like it. I say go for this version, but @jeromekelleher suggested gsl, so let's wait to see if he has something else to say.

[petrelharp]

Also, I also think that what you have proposed is correct, thus lowering the probability of error (but not to epsilon!).

[jeromekelleher]

Looks good to me! If this is simpler and we're sure it's correct then I say go for it. I guess we could consult TAOCP volume 2 as well, where I'm sure Knuth has some things to say about this problem (I don't have my copy on hand, though).

[jeromekelleher]

Actually, come to think of it, here's a version I wrote a while back which I think I based directly on Knuth.

[castedo]

By many reasonable and practical definitions of "sure" and "correct", I am sure the non-GSL function above is is correct.

Regarding the function you link to, I think it's a nice benefit to have a function like that where the random number gets passed in. That makes it conducive to having some unit tests that also act as documentation of corner cases. So I'll proceed with the non-GSL PR, use the more-correct code I mention above, but refactor to have a function interface more like what you provided. And I plan to add a few simple unit test, partly for documentation reasons.

However, on a critical note for the function you linked, given a simple case like ten 0.1 doubles, that function you link to will read past the end of the probability array, probably trigger an read access violation, if not, it will return an invalid index (greater than or equal to the size of the array). Here's a Juypter notebook in C to show this simple bug case:

https://nbviewer.jupyter.org/url/ref.castedo.com/temp/round_too_short.ipynb

I'm happy to explain more why the function fails in this case. At the end I've added some printf("%a", x) statements which I find very useful in understanding floating point calcs. That along with thinking of the numbers in base 2 or 16. I wish I could say I find http://ref.castedo.com/binspeak/ useful too, but really it's just a bit of silliness. :)

I did find a pdf of TAOCP v2 and looked into it quickly. I think I found the book section most closely related to this topic. I found a topic that is close, but not into the details of this specific algo.

[jeromekelleher]

Sounds good, thanks @castedo, let's go ahead with your version.

[castedo]

Excellent plan! I'll update the new PR, close the old PR, and un-draft the new PR when I finish the various changes discusses above.