Non-Limber FKEM fails catastrophically

[damonge]

I've seen the current FKEM method fail catastrophically in some cases. For instance, the following code:

import numpy as np
import matplotlib.pyplot as plt
import pyccl as ccl

# Redshift distribution
z = np.linspace(0, 4.72, 60)
nz = z**2*np.exp(-0.5*((z-1.5)/0.7)**2)

# Bias
bz = 0.278*((1 + z)**2 - 6.565) + 2.393

plt.figure()
plt.plot(z, nz*bz)
plt.xlabel(r'$z$', fontsize=16)
plt.ylabel(r'$b(z)\,\phi(z)$', fontsize=16)

# Power spectra
ls = np.unique(np.geomspace(2, 2000, 128).astype(int)).astype(float)
cosmo = ccl.CosmologyVanillaLCDM()
tracer_gal = ccl.NumberCountsTracer(cosmo, has_rsd=False, dndz=(z, nz), bias=(z, bz))
cl_gg = ccl.angular_cl(cosmo, tracer_gal, tracer_gal, ls,
                       l_limber=-1, p_of_k_a=cosmo.get_nonlin_power(), p_of_k_a_lin=cosmo.get_linear_power())
cl_ggn = ccl.angular_cl(cosmo, tracer_gal, tracer_gal, ls,
                        l_limber=1000, p_of_k_a=cosmo.get_nonlin_power(), p_of_k_a_lin=cosmo.get_linear_power())

plt.figure()
plt.plot(ls, cl_gg, 'k-', label='Limber')
plt.plot(ls, cl_ggn, 'r--', label='Non-Limber')
plt.loglog()
plt.legend()
plt.xlabel(r'$\ell$', fontsize=16)
plt.ylabel(r'$C_\ell$', fontsize=16)

yields the result

I wonder if this is related to how wide the redshfit distribution is in this case. Perhaps some of the internal FKEM parameters should be made tweakable to deal with this?

@paulrogozenski @xfangcosmo , thoughts?

[paulrogozenski]

The source of this error comes from the log-spacing of the chi array. Since the original chi array goes to zero, I am assuming some other minimum value, hard-coded to be 1.0e-6, over which to make a log-spaced spline in chi. Since this does not accurately sample the entire range of chi, we find this oscillatory behavior. I think for large bins, such as this, the minimum chi should be increased to reflect the n(z) distribution. In this particular case, you can directly set chi_min here to be 1.0 and the behavior is as expected. I will try to come up with a solution that can do this on-the-fly

[paulrogozenski]

This temporary fix should give you the following result

[damonge]

Oh, interesting. Thanks a lot @paulrogozenski , that was very fast.

Just adding these two lines after generating the N(z) worked for me:

z[0] = 0.001
nz[0] = 0.0

An automatic fix for it would be good indeed. I wonder how CCL is dealing with weak lensing tracers and nonlimber now, since by definition the lensing kernel should go to z=0.