non-Limber with FKEM gets scaled by fkem_chi_min

[ntessore]

I am trying to compute the Cls for super thin tophat n(z) with angular_cl() and non-Limber using FKEM. Somehow, my spectra end up being scaled (but not exactly) by fkem_chi_min, which I do not understand. Am I doing something wrong here?

lmax = 50
t = pyccl.NumberCountsTracer(
    cosmo,
    dndz=(np.linspace(0.0, 0.0025, 100), np.ones(100)),
    bias=(np.linspace(0.0, 0.0025, 100), np.ones(100)),
    has_rsd=False,
)
for fkem_chi_min in 1e-5, 1e-6, 1e-7:
    cl = pyccl.angular_cl(
        cosmo,
        t,
        t,
        np.arange(lmax + 1),
        l_limber=lmax,
        fkem_chi_min=fkem_chi_min,
        fkem_Nchi=100_000,
    )
    plt.loglog(np.arange(lmax + 1), cl)
plt.show()

[iagolops]

Hi @ntessore

I'm almost sure that this is because of your top-hat distribution. I tested using a similar one, but with zeros in other ranges of redshift, and I found a consistent result for the values that you tested.

nz = np.zeros(100)
z = np.linspace(0,0.01,100)
nz[(z<0.003) & (z>0.0005)] = 1

The N(z):

[nikosarcevic]

@ntessore is this now solved? if not, can you pelase let us know. otherwise please close the issue.

[ntessore]

Hi @nikosarcevic, I wouldn't consider it solved. If you have n(z) which don't start with zero, something weird happens in the FKEM code. There is a workaround (pad your n(z) with zero) but the weird thing that happens should probably not happen.

[nikosarcevic]

thank you. I do have another implementaiton of fkme waiting so after that is reviewed and merged, I will test you issue with the new code. sounds good?

[ntessore]

Sounds very good, thanks

[nikosarcevic]

@ntessore ok I have solved it.
I am going to think about it a bit more before comitting but I do know how to make sure this doesnt happen.

ETA: this will be in the new commit as part of the PR #1262

[nikosarcevic]

hey @ntessore and @iagolops i have pushed a commit to the exisitng pr (see dcb8230). this should solve your problems hopefully. can you perhaps clone this branch and try it yourself? note the slight changes in the api for nonlimber (to be conssitent with standard ccl)

ETA: how I checked it now

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

from pyccl.nonlimber_fkem.core import nonlimber_fkem
from pyccl.nonlimber_fkem.legacy_fkem import legacy_nonlimber_fkem

lmax = 50
ells = np.arange(lmax + 1)

cosmo = ccl.CosmologyVanillaLCDM()

z  = np.linspace(0.0, 0.0025, 100)
nz = np.ones_like(z)

t = ccl.NumberCountsTracer(
    cosmo,
    dndz=(z, nz),
    bias=(z, np.ones_like(z)),
    has_rsd=False,
)

# FKEM configs to compare
fkem_chi_mins = (1e-4, 1e-5, 1e-6, 1e-7, 1e-8)
NCHI = 100_000

p_of_k_a  = "delta_matter:delta_matter"
pk_linear = "delta_matter:delta_matter"
limber_max_error = 1e-2

fig, (ax_top, ax_bot) = plt.subplots(
    2, 1, sharex=True, figsize=(6, 8)
)

colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]

new_curves = {} 
for chi_min, color in zip(fkem_chi_mins, colors):
    # legacy implementation
    ell_limber_old, cells_old, status_old = legacy_nonlimber_fkem(
        cosmo,
        t,
        t,
        p_of_k_a,
        ells,
        lmax,
        pk_linear=pk_linear,
        limber_max_error=limber_max_error,
        Nchi=NCHI,
        chi_min=chi_min,
    )

    # new implementation
    ell_limber_new, cells_new, status_new = nonlimber_fkem(
        cosmo=cosmo,
        tracer1=t,
        tracer2=t,
        p_of_k_a=p_of_k_a,
        ell_values=ells,
        ell_limber=lmax,
        pk_linear=pk_linear,
        limber_max_error=limber_max_error,
        n_chi=NCHI,
        chi_min=chi_min,
    )

    n = min(len(cells_old), len(cells_new))
    ells_eff = ells[:n]

    new_curves[chi_min] = (ells_eff, np.array(cells_new[:n]))

    ax_top.loglog(
        ells_eff,
        cells_old[:n],
        linestyle="--",
        linewidth=1.5,
        color=color,
        label=rf"legacy, $\chi_\mathrm{{min}}={chi_min:g}$",
    )

    ax_top.loglog(
        ells_eff,
        cells_new[:n],
        linestyle="-",
        linewidth=3.0,
        color=color,
        label=rf"new, $\chi_\mathrm{{min}}={chi_min:g}$",
    )

chi_ref = fkem_chi_mins[0]
ells_ref, cells_ref = new_curves[chi_ref]

for chi_min, color in zip(fkem_chi_mins, colors):
    ells_i, cells_i = new_curves[chi_min]
    n = min(len(ells_ref), len(ells_i))
    ells_eff = ells_ref[:n]

    ratio = cells_i[:n] / cells_ref[:n]

    ax_bot.semilogx(
        ells_eff[1:],
        ratio[1:],
        color=color,
        label=rf"$\chi_\mathrm{{min}}={chi_min:g}$",
    )

ax_bot.axhline(1.0, linewidth=1.0, color="k")
fs = 16
ax_top.set_ylabel(r"$C_\ell$", fontsize=fs)
ax_top.set_title(r"FKEM: legacy vs new for different $\chi_{\rm min}$")

ax_bot.set_xlabel(rf"multipole $\ell$", fontsize=fs)
ax_bot.set_ylabel(rf"new / new($\chi_{{\rm min}}={chi_ref:g}$)")

ax_top.set_title(r"FKEM: legacy vs new for different $\chi_{\rm min}$")
ax_bot.set_xlabel(r"multipole $\ell$", fontsize=fs)

ax_top.legend(fontsize=9, ncol=2)
ax_bot.legend(fontsize=9, ncol=2)

plt.tight_layout()
plt.show()