Address repetitive Hankel transform code in utils and covariance modules

[bclevine]

Separating this issue out from #17.

Goal

Currently, we have two different Hankel transforms implemented in the code. The first is an FFTLog method implemented in utils/hankel_transform_1d.py, which is used in the lens magnification bias modeling. The second is a matrix operator method implemented in covariance/projection/hankel_transform.py, which is used in the covariance modeling.

We should do the following:

1. Benchmark the methods against one another.
2. Explore ways to reduce redundancy, improve performance, and combine methods.
3. Consider moving all methods into a separate Hankel transform module.

[bclevine]

Benchmarking

In this comment I'll post some preliminary results I get from benchmarking.

Projected transform:

The FFTLog and matrix methods agree:

But the matrix method exhibits some ringing. (This power spectrum is $\gamma_{t,LSS}$ from the magnification bias). I also noticed that it's quite difficult to fine-tune the matrix generation parameters for very small scales, and it can take a very long time. A few tweaks to the existing code can fix both these issues but at the cost of some performance, and there may also be other parts of the covariance code that the tweaks affect (I didn't explore it very thoroughly).

Now what about the speed? The matrix method is close to constant in time, while the FFTLog method should be O(NlogN) for the number of data points. Depending on the number of $\ell$ modes we use to sample the power spectrum, a different method will win out. But we get good results for 1000 or even 500 points in $\ell$, so the FFTLog version is probably faster for generating the datavector.

That said, there is a significant benefit to using the matrix form: it is much more efficient for computing blocks of the covariance matrix, where we would otherwise have to iterate over the FFTLog function. I didn't even try to test the timing on this because it seems obvious to me that the matrix method will outperform the FFTLog.

Spherical transform:

The spherical projections do not agree. Shown below is the $j_0$ transform for the nonlinear power spectrum, computed via FFTLog, matrix method, and also CCL as an external benchmark. I'm pretty sure this is a bug in the matrix method that came from the original code -- see spherical_correlation() in Dani's original repository. A comment mentions that the Bessel function should be order n+0.5, but that doesn't seem to be implemented in the code. There is a j_f quantity that is computed but never used. Since this function is never actually used anywhere in the code, I suspect that it was just never fully implemented and tested since it wasn't needed. We therefore should probably stick with the FFTLog-based integration for the spherical transform.