Inconsistencies in Theoretical Predictions Arising from Integration and Interpolation Methods in CCL

[nikosarcevic]

Hello, awesome people!

I've encountered a significant challenge related to the integration stability in the Core Cosmology Library (CCL) which I believe warrants our attention. This issue primarily concerns the discrepancies in theoretical predictions due to varying resolutions of redshifts, particularly observable through different outcomes in Fisher analysis and potentially affecting MCMC 2-point function analyses.

Problem Description

In my recent work, I've observed that changing the resolution over which the redshifts (N(z)) are defined leads to noticeable differences in Fisher contour predictions. Specifically, the differences emerge when (N(z)) is defined over the range (0 < z < 4) with varying resolutions, as illustrated by the following scenarios:

• Resolution 1: 700 points
• Resolution 2: 800 points
• Resolution 3: 900 points
• Resolution 4: 1000 points

For each resolution, I generated redshift arrays using np.linspace(0, 4, resolution), and subsequently ran Fisher analyses (using FISK, my Fisher analysis pipeline). The resulting contours differed significantly across the resolutions, including changes in direction which suggests a deeper underlying issue.
I have solved it by migrating to a different derivative method (not using the finite difference method but using the "stem" method -- see this repository for explanation).
However, it is still crucial to determine what is causing the inconsistencies.

To showcase this problem, I am including the Fisher contours for LSST Y1 and Y10 from FISK, using different resolutions of the redshift range. Note that the benchmark (in solid black line) is the previous SRD result.
We can see that, depending on the chosen resolution of the redshift range, the resulting contours will greatly differ. This is a very serious problem because it is impossible to gauge which result is the "true" one, even in case we have a benchmark to compare to.

Suspected Cause

My investigations lead me to believe that the issue may stem from a combination of numerical integration and interpolation strategies used within CCL, particularly those involving GSL integration methods and fallback routines. When certain methods fail or yield unstable results, CCL seems to switch to alternative strategies. While this ensures that a result is always produced, it might lead to unphysical outcomes due to inconsistencies in the integration or interpolation approach adopted. Integration, interpolation, and normalisation is used throughout CCL and I believe it is important to inspect and investigate the cause.

Impact

This problem is not just limited to Fisher analysis but extends to any computation relying on consistent theoretical predictions from CCL, such as MCMC analyses using 2-point functions, or any future anlyses from DESC. Inconsistent integration or interpolation methods could compromise the reliability of our results, making it imperative to address this issue.

Proposed Solution

To mitigate these discrepancies and ensure the robustness of CCL's theoretical predictions, I propose a thorough review and potential overhaul of the integration and interpolation methods used within the library. This could involve:

• A deep dive into the integration routines, especially those involving GSL, to identify potential sources of instability.
• Examination of the interpolation methods and their fallback strategies to ensure consistency across different resolutions.
• Development of a standardized approach or a set of guidelines for handling numerical integration and interpolation within CCL, aimed at minimizing inconsistencies.

I believe addressing this issue is crucial for the integrity of research utilizing CCL. I welcome any insights, suggestions, or collaborative efforts from the community to tackle this challenge together.

I look forward to your feedback and any further discussions on how we can improve CCL's reliability and accuracy.

If you encounter(ed) similar problems, please let me know and we can work on this together, write a report and solve this problem.

Best,
Niko

P.S. Let me know if there are typos or untrue statements in my post.

[ztq1996]

Hi there -

I have recently encountered a similar issue. I am using CCL to generate 3x2pt data vectors and test the sensitivity of cosmology to different photo-z parameters. The issue that threw me off is when my n(z) distribution is shifted, the cosmological bias jumps up and down at some threshold, see figure for an example below:

After some brief investigation, I realized the problem only occurs when I am adding outliers to my n(z). The problem seems to go away when I am using a simple Gaussian n(z), though I have no idea where is the culprit.

Investigating this would be really helpful, thank you!! Also happy to provide more details about my case if needed.

TQ

[nikosarcevic]

Hi there -

I have recently encountered a similar issue. I am using CCL to generate 3x2pt data vectors and test the sensitivity of cosmology to different photo-z parameters. The issue that threw me off is when my n(z) distribution is shifted, the cosmological bias jumps up and down at some threshold, see figure for an example below:

After some brief investigation, I realized the problem only occurs when I am adding outliers to my n(z). The problem seems to go away when I am using a simple Gaussian n(z), though I have no idea where is the culprit.

Investigating this would be really helpful, thank you!! Also happy to provide more details about my case if needed.

TQ

Many thanks for reporting this!

I think other folk have more examples so it would be great if we can collectively figure this out and solve it.

[nikosarcevic]

@ztq1996 so basically your problem goes away after you smooth things out? it definitly smells like the problem we think it is causing this -- the interpolation in ccl.

[carlosggarcia]

Sharing here what I commented in today's MCP telecon: could it be that by increasing the binning of your N(z) by a lot, they become noisy? (This is just a guess based on @nikosarcevic presentation)

[nikosarcevic]

Many thanks for sharing ideas.
I personally think it is not the case as the accuracy should increase. This is not resolution of 10000000, it is just 800 or 100 so not much.
I wanted to make more plots, for example, the impact onto the data vector, that might be interesting to see.

The fact is that other people are seeing this as well (at various stages and various outputs), like TQ above and also we became aware of it during the CCL testing/benchmarking. Marianna PL and her student reported it (I will ask them to share their finding here). I think we can also ask @marcpaterno @vitenti and Javi to chip in here as they also said they noticed this problem.

[fjaviersanchez]

If you increase the binning of a real distribution it does become noisy, but I think that @nikosarcevic is using an analytical form. for the dndz (please correct me if I am wrong!)

I think this could be an issue of the spline interpolation being noisy for certain choices of binning, but at the same time, this opens the question of what we are going to do with the data. As you mention @carlosggarcia the binning there becomes pretty important, but I guess that might be a question for the PZ WG :).

[nikosarcevic]

yes i forgot to answer here -- indeed im using an analytical form

[damonge]

Have you plotted the C_ells themselves? My guess is that there is some numerical instability there, and this is what causes the issue. It'd be good to look at those first to find a solution.

[nikosarcevic]

yeah i didnt have the time but i mentioned yesterday I will plot the data vectors and share here

[nikosarcevic]

Hi all, @fjaviersanchez @damonge et al

sorry for waiting, I had to eat

I quickly made a C_ell comparison for different z resolutions. I only checked one cross and one auto correlations: tomo bins (0,1) and (0,0)

Here are the results: we see the spikes. it is small but maybe significant overall?

[damonge]

Yes, those look terrible. It'd be good if you can plot:

• The N(z)s you're passing in each case.
• The corresponding lensing kernels (you can get those from the tracers with get_kernel().
  These are the two things that determine the C_ell

[nikosarcevic]

Yes, those look terrible. It'd be good if you can plot:

• The N(z)s you're passing in each case.
• The corresponding lensing kernels (you can get those from the tracers with get_kernel().
  These are the two things that determine the C_ell

haha, on it already! working on N(z) atm, will add kernel as well
many thanks, David!!

[nikosarcevic]

@damonge et al. here it is (sorry, i couldnt do it sooner)

getting kernels from sources (tracers) in bin 1 and 2 (low) and bins 4 and 5 (high z).
Like in previous plots for the data vectors, I divide each W with the first one. Each W is obtained for a different z resolution. So I divide W(res=300)/W(res=300), W(res=300)/W(res=500), W(res=300)/W(res=800), etc.

seems that things are much worse at low z which is a concern.

N(z) is fine (comes from my code, modified version of it is on cclx)

also, @paulrogozenski has been kind enough and joined the effort. he is also investigating this and will report (I think)

[damonge]

OK, thanks. I'm not sure I understand what's in the X axis of the W plots, and also unsure what normalisation you're using. Can you just plot W(chi), as it comes out of the get_kernel function, without any normalisation, for different resolutions?

[nikosarcevic]

Hi, David

ok yes it is very confusing and I am not really sure I completely understand it. But this is what I know:

Here is the snippet of how I obtain kernels W:

# Create WeakLensingTracer objects for the two source bins
        source_1 = ccl.WeakLensingTracer(cosmology,
                                         dndz=(self.redshift_range, source_bins[key_1]),
                                         ia_bias=(self.redshift_range, ia_source))
        source_2 = ccl.WeakLensingTracer(cosmology,
                                         dndz=(self.redshift_range, source_bins[key_2]),
                                         ia_bias=(self.redshift_range, ia_source))

        #chi = ccl.comoving_radial_distance(cosmology, 1 / (1 + self.redshift_range))
        chi = None
        kernel_1, chi_values_1 = source_1.get_kernel(chi=chi)
        kernel_2, chi_values_2 = source_2.get_kernel(chi=chi)

so, calling get_kernel() on each tracer will give me a tuple.
If chi=None, each item in the tuple will be a list of 2 arrays.
In other words, kernel_1 will be a list with two arrays aka kernel_1 = [(kernel_1_array_1), (kernel_1_array_2)]. Same for chi.
Shape wise (when I plot kernel_1_array_1 and kernel_1_array_2), they look the same, only kernel_1_array_2 seems to be normalized (they differ in scaling/amplitude).

In case I pass chi = ccl.comoving_radial_distance(cosmology, 1 / (1 + self.redshift_range)) instead of None, I do not get a list of two arrays for tuple items, but only this kernel_1 = [(kernel_1_array_1)].

If I look at the docstring for get_kernel(), this is what I see

Help on function get_kernel in module pyccl.tracers:

get_kernel(self, chi=None)
    Get the radial kernels for all tracers contained
    in this ``Tracer``.
    
    Args:
        chi ((:obj:`float` or `array`)): values of the comoving
            radial distance in increasing order and in Mpc. If ``None``,
            returns the kernel at the internal spline nodes.
    
    Returns:
        `array`: list of radial kernels for each tracer. The shape
        will be ``(n_tracer, chi.size)``, where ``n_tracer`` is the
        number of tracers. The last dimension will be squeezed if the
        input is a scalar.

[BTW, I find it confusing what is meant by "The last dimension will be squeezed if the
input is a scalar." the last sentence in the docstring.]

So, plotting again would be redundant as I already have it. Judging from the docstring, I assume that x axis is n_tracer?
I do not know if this helps. At one point last night, Paul jumped in with me on Zoom to try to figure out but we didnt really get anywhere.

Does this help at all?

[damonge]

Hmm, maybe I can play around with this. Do you mind giving me the analytical dndz you're sampling here?

[nikosarcevic]

Most certainly!

The code from my pipeline to generate the SRD N(z) and tomobins is actually available on CCLX:

• SRD N(z)
• SRD tomo bins

[nikosarcevic]

Also if it is of any help, here is how I sample over different z resolutions. This is calling my pipeline but maybe it is better you can have a clear picture of how I did this

ETA: it is not the best code! apologies! but hope it helps somewhat

%%time
# My redshift ranges
resolution_list = [300, 500, 800, 1000, 1300, 1500, 2000]
redshift_ranges = [
    np.linspace(0.01, 4., resolution_list[0]),
    np.linspace(0.01, 4., resolution_list[1]),
    np.linspace(0.01, 4., resolution_list[2]),
    np.linspace(0.01, 4., resolution_list[3]),
    np.linspace(0.01, 4., resolution_list[4]),
    np.linspace(0.01, 4., resolution_list[5]),
    np.linspace(0.01, 4., resolution_list[6]),
]

# Creating presets list
presets_list = [
    Presets(analysis="cosmic_shear", forecast="srd", forecast_year="1", redshift_range=redshift_range)
    for redshift_range in redshift_ranges
]

# Defining cross and auto pairs upfront
cross_pairs = [(3, 4)]
auto_pairs = [(4, 4)]

# Initializing a list to hold results for each preset
results_list = []

for presets in presets_list:
    # For each preset, calculate source bins and ells
    sources = RedshiftDistributionModel(presets).source_sample_redshift_distribution()
    source_bins = Binning(presets).source_bins()
    ells = MultipoleScale(presets).get_ells()
    
    cs = CosmicShear(presets)
    
    # Initialize dictionaries to hold cls values and kernels for cross and auto pairs
    cls_cross = {}
    cls_auto = {}
    
    kernels_cross = {}
    kernels_auto = {}
    
    # Calculate cls and kernels for cross pairs
    for i, j in cross_pairs:
        cls_key = f"{i}_{j}"
        cls_cross = cs.cosmic_shear_cls(i, j)
        kernel_1_cross, chi_1_cross, kernel_2_cross, chi_2_cross = cs.get_kernels(i, j)
        kernels_cross[f"kernel_1"] = kernel_1_cross
        kernels_cross[f"kernel_2"] = kernel_2_cross
        kernels_cross[f"chi_1"] = chi_1_cross
        kernels_cross[f"chi_2"] = chi_2_cross
    
    # Calculate cls and kernels for auto pairs
    for i, j in auto_pairs:
        cls_key = f"{i}_{j}"
        cls_auto = cs.cosmic_shear_cls(i, j)
        kernel_1_auto, chi_1_auto, kernel_2_auto, chi_2_auto = cs.get_kernels(i, j)
        kernels_auto[f"kernel_1"] = kernel_1_auto
        kernels_auto[f"kernel_2"] = kernel_2_auto
        kernels_auto[f"chi_1"] = chi_1_auto
        kernels_auto[f"chi_2"] = chi_2_auto
    
    # Convert cls to numpy arrays for consistency
    cls_cross = np.array(cls_cross).T
    cls_auto = np.array(cls_auto).T
    
    # Store the results in a structured way
    results_list.append({
        "cls_cross": cls_cross,
        "cls_auto": cls_auto,
        "ells": ells,
        "sources": sources,
        "source_bins": source_bins,
        "kernels_cross": kernels_cross,
        "kernels_auto": kernels_auto,
    })

[damonge]

I can't seem to reproduce the issue. This is the result I get

with the following code:

cols = ['k', 'r', 'g', 'b', 'y', 'c', 'm']
nsamples = [300, 500, 800, 1000, 1300, 1500, 2000]

cosmo = ccl.CosmologyVanillaLCDM()
ls = np.unique(np.geomspace(2, 2000, 128).astype(int)).astype(float)
dndzs = []
cls = []

for nsamp, col in zip(nsamples, cols):
    z = np.linspace(0.01, 3.5, nsamp)
    dndz = SRDRedshiftDistributions(z, 
                                    galaxy_sample="source_sample",
                                    forecast_year='1').get_redshift_distribution(normalised=True)
    dndzs.append((z, dndz))
    t = ccl.WeakLensingTracer(cosmo, dndz=(z, dndz))
    cl = ccl.angular_cl(cosmo, t, t, ls)
    cls.append(cl)

plt.figure()
for nsamp, cl, col in zip(nsamples, cls, cols):
    plt.plot(ls, cl/cls[0]-1, col+'-', label=f'{nsamp} samples')
plt.legend()
plt.xscale('log')
plt.xlabel(r'$\ell$', fontsize=15)
plt.ylabel(r'$\Delta C_\ell/C_\ell^0$', fontsize=15)

The differences are at the level of 1E-4, which is the level of accuracy that CCL reaches. If differences of this kind actually impacted forecasts so much, we'd be in huge trouble (since astrophysical uncertainties are certainly much, much larger than that!).

My guess is that the issue is somewhere else (e.g. accuracy of numerical derivatives, if this is Fisher-based? Has anyone tried running a quick MCMC?).

[nikosarcevic]

The initial inconsistencies are in Fisher, yes. I have solved the fluctuation problem by implementing a different derivative method (I linked it in the original post above).

However, other people are reporting instabilites in non-fisher related stuff. I will ask Marianna PL to explain from her side.
Paul will share his version of the code, as well as an old issue he raised on ccl slack re gsl.
Also, I will as TQ to exaplain further the problems he is seeing.

No one tried MCMC at this point but we could and should do it. (I cannot do it now, I am also working on the baryons halo model with Elisa)

@fjaviersanchez @vitenti -- can you two let us know what you saw regarding this issue? do you have any plots?

[nikosarcevic]

@ztq1996 woudl you be so kind and maybe give us a bit more clarification and description of your problem? ;grateful:

[luigilcsilva]

Hi, everyone!

Some time ago, I was computing some auto angular power spectra and cross angular power spectra for galaxies, using the bin1_histo.txt and bin2_histo.txt dndz's, which can be found here. I was comparing the results from CCL and NumCosmo. In the case of CCL2 cosmology, using the bin1_histo distribution and computing the auto angular power spectrum, I end up with some problematic points. For l = [ 519 586 730 798 931 1063 1147 1148 1189 1190], the accuracy value was above the expected limit of 0.1. In a further investigation, I found that the Cls as computed by NumCosmo are actually smooth in these regions, while the Cls as computed by CCL have some strange fluctuations. We also think this may have something to do with the interpolation and integration methods used in CCL.

[luigilcsilva]

Hello, Niko!

Yes, we think that the kinks in the N(z) and the way that CCL is dealing with them when interpolating may be impacting in the Cl computed values. One of the next steps we want to do is to plot the CCL and NumCosmo interpolated dndzs directly, and compare them.

[nikosarcevic]

many thanks!

may I also ask if you would be so kind to check for example if you pass a smooth N(z) into ccl and see what the resulting Cl looks like?

[luigilcsilva]

You are welcome!

Well, I did this for gaussians generated from the same bin1_histo and bin2_histo dndz's:

For the auto angular power spectra, using the gaussian generated from the bin1_histo N(z), and for CCL2 cosmology, the accuracy has a very good behaviour. In the image below, there are the accuracy plots, using n=[100, 200, 500, 1000] number of bins in the gaussian dndz distribution.

Because the accuracy had a good behaviour, I didn't zoomed the Cl plot as I did in the previous case.

The accuracy has a good behaviour for CCL1 and CCL3 cosmologies as well (I tested just these three, CCL1, CCL2 and CCL3).

[nikosarcevic]

great!
ok, i need a bit more help:

1. what is the A_acc -- would you be so kind and explain what is this quantitiy and how you calculate it?
2. what is the orange line?
3. why is the blue result resonant?

(forgive my ignorance im just trying to understand )

[luigilcsilva]

No problem! Sorry I didn't explain it earlier.

1. The A_acc is the accuracy metric for the comparison. It is given by equation 95 and 96 of the CCL 2019 article (page 28). I use C^{CCL}_l computed from CCL, and I use the values from NumCosmo for the reference parameters Cl and \sigma_l.

2. The orange line is the accuracy requirement, also according to the CCL 2019 article, just below equation 96 (page 28). Quoting the article, "differences must be smaller than one tenth of the cosmic-variance errors", for the Cl's.

3. For this question we still don't have an answer. In the gaussian dndz's case, these oscilations are an order of magnitude below the accuracy requirement, so we are less worried. But yes, it's worth investigating this. Looking at the accuracy metric, it probably has to do with the differences in the computation of the Cls in both libraries.

[fjaviersanchez]

Apologies for the slow response. I do see differences, but in my case it's because the way I am generating the dndz itself depends on the the sampling that I do in dndz. The first plot is the relative difference in C_ells and the second plot it's the absolute difference in the dndz.

And the code:

from augur.tracers.two_point import SourceSRD2018
cols = ['k', 'r', 'g', 'b', 'y', 'c', 'm']
nsamples = [300, 500, 800, 1000, 1300, 1500, 2000, 10000]

cosmo = ccl.CosmologyVanillaLCDM()
ls = np.unique(np.geomspace(2, 20000, 1000).astype(int)).astype(float)
dndzs = []
cls = []

for nsamp, col in zip(nsamples, cols):
    z = np.linspace(0.01, 3.5, nsamp)
    srcs = SourceSRD2018(z, Nz_nbins=5, Nz_ibin=0, Nz_alpha=0.78, Nz_z0=0.13, Nz_sigmaz=0.05)
    dndz = srcs.Nz
    dndzs.append((z, dndz))
    t = ccl.WeakLensingTracer(cosmo, dndz=(z, dndz))
    cl = ccl.angular_cl(cosmo, t, t, ls)
    cls.append(cl)

plt.figure()
for nsamp, cl, col in zip(nsamples, cls, cols):
    plt.plot(ls, cl/cls[0]-1, col+'-', label=f'{nsamp} samples')
plt.legend()
plt.xscale('log')
plt.xlabel(r'$\ell$', fontsize=15)
plt.ylabel(r'$\Delta C_\ell/C_\ell^0$', fontsize=15)

plt.figure()
from scipy.interpolate import interp1d
spl = interp1d(dndzs[0][0], dndzs[0][1], kind='cubic')
for i in range(len(dndzs)):
    plt.plot(dndzs[i][0], dndzs[i][1]-spl(dndzs[i][0]), label=f'{nsamples[i]} samples', c=cols[i])
plt.legend()
plt.xlabel('$z$')
plt.ylabel(r'$\Delta dndz$', fontsize=15)

I'll try with some other dndz more numerically stable :)

[fjaviersanchez]

This is coming from the augur side because I am sampling the overall dndz, and then I integrate and convolve it depending on the number of bins that I choose.

[fjaviersanchez]

I checked with another fun dndz that should be easy to deal with for CCL:

dndz = (np.exp(-(z-0.5)**2/(0.05*(1+0.5))**2) + np.exp(-(z-0.45)**2/(0.03*(1+0.45)**2))*0.05
            + np.exp(-(z-0.7)**2/(0.1*(1+0.45)**2))*0.01)

Same code as above just changing dndz for that. And I do see differences in the C_ells but of order 10^-5:

The C_ells are more different as the dndzs at low-z are more different, which makes sense. I think these are just very hard to avoid numerical issues due to having to numerically integrate the power spectrum and the dndz.

[nikosarcevic]

@fjaviersanchez so you think this is fine then?

[fjaviersanchez]

I think that for Fisher this is probably fine(?) as we know that it has certain limitations. We can always say: "we chose this prescription" or we know that with this prescription the results are robust to x% on certain parameters. For other things, I'd say it is also probably fine because there will be effects with a larger impact (I think that for certain scales the models are much more uncertain that this), but we'll need to evaluate in a case-by-case basis, I guess.

[nikosarcevic]

yeah i also feel we should check if it propagates and if yes how it propagates to other objects

i tihnk waht @paulrogozenski raised is interesting and we do not want to be really bad in low z. at least that is my impression

[paulrogozenski]

I'm finding that for multiple bins, the resolutions of your dndz sampling matters more for low-z bins than for high-z bins. In the low-z regime, a low number of z samples in your dndz distribution changes the bounds of limber integration relative to a higher-resolution sampling. These differences are more catastrophic at low-z, as this changes the comoving distance more than at high-z. I've attached plots showing the ratios of Cls between low and high resolution dndz for low-z (lens0) and high-z (len4) autocorrelation. Chi_min and chi_max for lens0 can change as much as 5% each between the low-resolution and high-resolution dndzs, but less than 1.5% for lens4.

[nikosarcevic]

@paulrogozenski you knwo what im thinking? cna you do the same plots but with 255, 256 and 257 samples/resolution?
there is this 256 limit if i remember

[paulrogozenski]

Sorry I missed this.The 256 sampling is not a limit, but rather the default resolution of the WeakLensingTracers and NumberCountsTracers. Varying around this default value, though, I find that the lens bins for these schemes are consistent with one another, but the low-z source bins are not.