Faster prior variance calculation

[nspope]

Use a recursive approach to calculating prior variances that is a couple of orders of magnitude faster:

so it's possible to get the exact variances for 100,000 tips in a few minutes (by comparison, it's infeasible to get exact variances for more than 2,000 tips or so with the original algorithm).

Some indication that the two algorithms are equivalent, showing sqrt(variance) for 1280 tips:

Still trying to get my head around the interpolation scheme used for the approximation, and if this algorithm will be useful there (it recurses backwards from the maximum number of tips).

[nspope]

Some interesting (but possibly useless) applications of this algorithm:

• Get logarithmic moments -- this would allow matching lognormal/gamma via sufficient statistics rather than mean/variance
• Adjust for variable population size at the level of the hypoexponential components (e.g. do the time rescaling before marginalizing over number of extant ancestors, and mixing over node spans) -- this should give more accurate moments for the gamma/lognormal approximations, but maybe not enough to prioritize?

[hyanwong]

Really neat. Thanks for this Nate.

[nspope]

Darn, this is working but fried some unit tests.

[nspope]

This is good to go. I moved the original, slow variance calculation into utility_functions and test vs fast implementation.

[hyanwong]

LGTM. Maybe it's worth keeping a few tests for the simpler version? Also, using the approximation only when N>20,000 seems quite a high threshold. I assume it's so fast that this is fine? And I assume that's when you removed the "Calculating Prior Variances" tqdm feedback?

[hyanwong]

Also, for very large inferences, if e.g. N=1e6 and the user has set approx_priors=False, is it worth outputting some sort of progress bar? Perhaps we only display a progress bar for this stage when progress=True and approx_priors=False?

[nspope]

Hmm, I removed the progress bar because the variance calculation is done recursively in numba -- so, no way to output progress without slowing things down dramatically.

[nspope]

using the approximation only when N>20,000 seems quite a high threshold

I think it takes around 1/2 a minute, according to the benchmark above? That's really when the cost starts kicking in.

[hyanwong]

Hmm, I removed the progress bar because the variance calculation is done recursively in numba -- so, no way to output progress without slowing things down dramatically.

Oh, well, if there's no obvious trick there then just drop it.

We could output a logging warning that this might take a while if approx_priors=False and N is so large that it will take of the order of an hour or more?

[nspope]

Yes, that's a good idea

[nspope]

Ok I think this is done-- except I'm not sure about the default for using the approximate prior. I suppose that if there's missing data, the calculation has to repeated for every realized # of tips across the tree sequence? So, I set the default to 10,000 which should take on the order of a second according to benchmark above. Seem OK?

[hyanwong]

A second sounds OK to me, as does a default of 10,000, so I'll just merge this. Thanks a lot Nate!