ExponentialGrowth numerical issues

[tgvaughan]

The current implementations of the getIntensity() and getInverseIntensity() methods for the exponential growth population model are numerically unstable for extremely small but non-zero growth rates.

For instance, here is the result of computing getIntensity() for times between 0 and 10, a present day population size of 1 and three different growth rates. All of these curves should be indistinguishable from one another:

It's easy to see why this problem occurs by looking at the definition of getIntensity:

    public double getIntensity(double t) {
        double r = getGrowthRate();
        if (r == 0.0) {
            return t / getN0();
        } else {
            return (Math.exp(t * r) - 1.0) / getN0() / r;
        }
    }

While the limit converges to t/N0 as r goes to 0, when r is small but finite the code above begins dividing a very small number by another very small number.

While this instability probably rarely causes issues, it can (does) introduce difficulties when trying to simulate under such distributions. As coalescent simulations sometimes occur in operators, instabilities like this can cause bad things to start happening to the sampler.

Ideally we would come up with a more numerically stable implementation of this function and the associated inverse intensity. Another could be to directly implement a stable method to handle the common use case for these functions, i.e. computing getInverseIntensity(getIntensity(t) + delta) for some delta.

[rbouckaert]

Ouch -- I like the getInverseIntensity solution for this.