LSSTDESC / CCL

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

support for 3d xi(r) and wp(rp) #219

Open jablazek opened 7 years ago

jablazek commented 7 years ago

CCL should be able to calculate and output 3d correlation function xi(r) and the projected (physical) w_p(r_p). xi(r) is very straightforward with P(k) and FFTLOG. w_p(r_p) is similarly simple in the Limber approx.

tmcclintock commented 7 years ago

@jablazek can you link the repo for FFTLog? I am an idiot and haven't located it yet.

tmcclintock commented 7 years ago

Just something to note, but another fast algorithm exists to do this integral using a quadrature rule. I have implemented it here, based on this algorithm.

jablazek commented 7 years ago

i assume that they were using @slosar version: https://github.com/slosar/FFTLog

tmcclintock commented 7 years ago

Now that fftlog.c is in the master branch at src/fftlog.c this amounts to just linking and making an interface somewhere, or else letting that file be the interface itself. @jablazek thoughts?

jablazek commented 7 years ago

@elisachisari @tmcclintock @sukhdeep1989

Agreed, modulo the modifications made to fftlog to operate on projected statistics.

sukhdeep2 commented 7 years ago

Computing w_p is in principle straight forward. It is the same transform as going from Cl->w(\theta). Right now, the ell limits are hard coded which is causing an issue to transform P(k) since k range is very different.

sukhdeep2 commented 7 years ago

3D xi(r) will require separate function since the bessel function is different. This is actually easier since FFTlog is coded to compute xi(r) and we had to hack around it.

jablazek commented 7 years ago

w_p is super easy in the limber approx.

sukhdeep2 commented 7 years ago

These are now implemented in correlation branch. Code compiles (C part) but output needs to be tested.

elisachisari commented 6 years ago

xi(r) was addressed by #306 and merged.

beckermr commented 5 years ago

Note for future ppl. Only wp(rp) remains on this issue.

jablazek commented 3 years ago

Sprinting on this with @Andromedanita. Will try to finish implementation tomorrow. We need a new function similar to ccl_correlation_3d in ccl_correlation.c.

The only difference should be the first argument in ccl_fftlog_ComputeXi3D. In the standard 3d correlation call, l = 0, corresponding to mu = l +0.5 for a J_mu transform.

We may also need to apply the correct prefactor. Also how are k and r corrections applied in this fftlog implementation? (what we call the alpha and beta factors)

Pinging @sukhdeep1989 .

sukhdeep2 commented 3 years ago

@jablazek I didn't fully understand the question. You mean prefactor for k,rp, similar to those of ell, theta?

For w(rp) like things, here's another proposal: Why don't we implement the equations from Baldauf's paper that we once talked about. I wrote a generalized form of those eqs. here. (eq. 14-18). These include the RSD corrections as well and then we may not always need non-limber ;), since RSD is the dominant correction. For wg+ like terms, there is an additional step to convert wg_kappa to wg+.

jablazek commented 3 years ago

@sukhdeep1989 : my question about the prefactor and the k, r scalings was in order to allow the generic bessel tranform using fft-log. (e.g. need additional k factors to convert between integer and half-integer bessel functions). I'm sure this is somewhere in the fft-log implementation.

If I understand those equations correctly, it is slightly orthogonal issue. We both want to allow for an anisotropic Pk and to perform e.g. the J2 integral to get wp(rp). But you are saying that these two things should be done at the same time to allow for non-limber?

sukhdeep2 commented 3 years ago

@jablazek I had a quick looks at ccl_fftlog_ComputeXi2D and ccl_fftlog_ComputeXi3D in https://github.com/LSSTDESC/CCL/blob/master/src/ccl_fftlog.c ? As far as I can tell, no prefactors are needed. The two functions pass a different number for dimension (2D vs 3D) and ccl_fftlog_ComputeXi3D adds 0.5 to the order of bessel function. They get the respective correlation functions back in return. There are additional beta dependent corrections applied to the correlation function multipoles in ccl_correlation_multipole, inside ccl_correlation.c .

On the second point, what I'm trying to say is that we can do wp this way to apply RSD like anisotropic corrections. Remember IA also has such terms (1-mu^2 => beta=-1 in RSD notation). Doing a direct J0, J2 like integrals miss them. It is not full non-limber, it only includes RSD which is the dominant correction (and it is flat sky). To do a J0 equivalent, one can simply pass beta=0.

Implementation wise, we can follow the existing xi_multipole functions (ccl_correlation_multipole). wp will require a different set of prefactors on multipoles than what is used inside ccl_correlation_multipole. I guess we can do both types of implementations, since simple J0,J2 should be easy with ccl_fftlog_ComputeXi2D .