LSSTDESC / CCL

DESC Core Cosmology Library: cosmology routines with validated numerical accuracy
BSD 3-Clause "New" or "Revised" License
133 stars 64 forks source link

correlation_3d for fixed scale factor instead of using Pk2D #1176

Open James11222 opened 2 months ago

James11222 commented 2 months ago

Hello,

I've been using the halos module for my work and making some modifications for my own research cases. Currently I have a question pertaining to computing correlation functions with the halos module.

In the halos.pk_2pt module it is fairly easy to compute the power spectrum for a fixed scale factor using halomod_power_spectrum() allowing me to avoid creating a Pk2D object which is done with halomod_Pk2D(). I would like to be able to compute the 3D 2pcf ξ(r) from this power spectrum at a fixed scale factor without having to compute the power spectrum for several other scale factors, but when I go to use the correlations module and use correlations_3d(), it only takes a Pk2D object as an argument for power spectra. I've tried making a Pk2D object with a single scale factor but I get that the kernel crashes in this case.

Below I include an example snippet of what I'm talking about, and what I would like to do

#compute pk at a fixed scale factor using the halo model
P_mm_ccl = ccl.halos.halomod_power_spectrum(cosmo, hmc, k=k_arr, a=np.array([a_sf]), prof=prof_NFW) 

#convert pk to a Pk2D object with fixed scale factor
P_mm_ccl_2D = ccl.pk2d.Pk2D(a_arr = np.array([a_sf]), lk_arr = np.log(k_arr), pk_arr=P_mm_ccl) 

#compute the 3D 2pcf from the Pk2D object
xi_mm_3D = ccl.correlation_3d(cosmo, r=r_arr, a=a_sf, p_of_k_a=P_mm_ccl_2D)

but this crashes the kernel regardless of machine I'm using. Ideally I'd like to bypass the second step altogether since my analysis is all done at a single redshift.

Is there a way to compute the 3D 2pcf from a power spectrum at a fixed scale factor rather than with a Pk2D object requiring multiple scale factors? To my understanding this doesn't seem to be implemented at the moment.

Onwards, James

damonge commented 2 months ago

@James11222 indeed this is not implemented. The easiest short-term, "hacky" solution for you would be to create a Pk2D that spans more than one scale factor. Since you don't care about evolution, you can probably just repeat the same Pk array at a few values of a (I forget what kind of interpolation we're using here, but I imagine 3 values of a will be enough). Let me know if this works for you.

The longer-term solution would be to allow for Pk2Ds to be generated at a single scale factor (and assume that they are the same at all as) or, perhaps even simpler, to allow for factorisable Pk2Ds, so you can pass the k and a dependence as two different multiplicative factors (this is already implemented at the C level, but not visible to python). If this would be a desirable feature, perhaps best to open a dedicated issue for it.

tilmantroester commented 2 months ago

We should catch the C-level error in any case though, crashing the python interpreter is bad.