Possible error in implementation of FFTLog algorithm #2
Description
Hi hankl-team,
recently I have worked on my own implementation of the FFTLog algorithm and I compared my results with your code (in particular the function P2xi()). I came to notice that there are tiny differences. I consequently looked into the source code of hankl and found a possible error in the _u_m_term() function. Here you define omega = 1j * 2 * np.pi * m / L as is done in the original paper by Hamilton (https://arxiv.org/abs/astro-ph/9905191, Eq. 174). This expression uses the fact that the input function is sampled on ln(r_n) =ln(r_0) + n*Delta_L and therefore 1 / (N*Delta_L) = N / (N*L) = 1/L due to the fact that Delta_L = L/N. L is hereby the logarithmic period and N the number of sampling points.
However, in that article a summation from -N/2 to N/2 is assumed in the discrete Fourier transform (DFT) which would result in a total of N+1 sampling points. Hence, the N/2 and -N/2 contributions are weighted by 1/2. In your implementation the DFT is a FFT that sums only from 0 to N-1 and hence uses N sampling points. This in turns implies that Delta_L = L / (N-1) and therefore the term for omega needs to be corrected to omgea = 1j * 2 * np.pi * m / L * (1. - 1./N). I changed the hankl code accordingly and checked both the analytic example you provide in the docs and the cosmology example comparing the original hankl version with the fixed one but also to an external package for the cosmology example.
For the analytic example is used r = np.logspace(-7, 7, N) with N being in [256, 512, 1024], which is a bit different from the example in the docs. As you can see, the absolute difference between the analytic function and the output of the FFTLog is orders of magnitude smaller when the error is fixed. I cannot exactly explain why the difference in the fixed version for N=1024 is a bit larger than the N=512 case but if I enlarge L, at some point the difference is the smallest for N=1024.
I also tested the cosmology example from the docs where I set k = np.logspace(-6, 1.3, 2**10) and just transformed the linear real-space power spectrum into configuration space. I then interpolate the output on s = np.linspace(10, 200, 100) using scipy's interp1d function. I also compare the fixed version of hankl to mcfit (https://github.com/eelregit/mcfit) yielding a very good agreement. Note that the effect seems to be present on scales around the BAO possibly affecting the recovery of certain cosmological parameters.
I also found the same error in the definition of the low-ringing condition where you define Delta_L=L/N that should instead be a Delta_L=L/(N-1). It seems like fixing this bug yields again very good agreement with mcfit once lowring=True for both. Although, looking at the difference between the green and the blue curve in the lower panel, the impact of the error on the low-ringing condition is much smaller than on the _u_m_term()-function.
I am happy to provide the code to produce the plots and I am looking forward to hear back from you.
All the best,
Martin