fast P(k,z) access

[plaszczy]

Dear Class support,
I would like to gain on the speed of P(k,z) computations (which is too low by default) and wonder how to do that efficiently so I have a few questions:

• is it possible to specify externally the k node values used ?. Is there a dedicated function for that? (is it good idea? it may depend interpolation scheme used)
• is a call to Pk(vec k,z) faster than individual P(k_i,z) ones?
  -there is a z_max_pk (that I set) : is there a z_min_pk one?
  -how is performed the P(k,z) evolution in z? (in which part of the code?)

any other idea?
thanks, my best
stephane

[ThomasTram]

Hi Stephane

Are you calling CLASS through the Python wrapper, the C++ wrapper or directly through C?

Cheers,
Thomas

[plaszczy]

Hi Thomas,
throught the standard C API.
I try to understand how p(k,z) is computed; actually what I need is fast access to p(k,z) but for a fixed z . it looks to me that this call performes a 2D spline interpolation whicle I only need 1D (in k)
Is there some explanation somewhere on this module?
thanks
stephane

[ThomasTram]

Hi again

Computing P(k,z) requires most of the CLASS modules except for the transfer module (that computes transfer functions in harmonic space) the lensing module (that computes CMB lensing) and the non-linear module (unless halofit is required). This takes of the order 1 second. After this computation the P(k,z) can be recovered using spectra_pk_at_z(), but the time needed for doing spline interpolation in this function call is completely negligible compared to running CLASS. But if you want, you can directly access P(k,z) for the computed (k,z)-values using the array psp->ln_pk[].

A real speedup in this computation can be obtained by
*) Not asking for Cls, in order to skip the transfer module.
*) Reducing the number of computed k-values, especially at large k.

For the last part, you can play with the settings mentioned in this issue #32.

Cheers,
Thomas

[plaszczy]

thanks Thomas,
actually my limitation really comes from calling spectra_pk_at_k_and_z not from computing the table ( I have many many calls). Is there a way to speed up the interpolation?
cheers
stephane

[ThomasTram]

Hi again

Okay, you must have many k-values then! But then spectra_pk_at_k_and_z() is indeed pretty slow, since for each k, the k-array at the given tau is constructed from interpolation and then splined. You will need to write a new fast function, but first a question: Do you know all the k-values you want to compute P(k) from the beginning, or do you generate them on the fly? (The latter would be the case if you are for instance doing adaptive integration in k...)

Cheers,
Thomas

[plaszczy]

right, this is for tomography where up to 10^10 P(k,z) values may be required. But good news k and z
are fixed and known
cheers
stephane

[ThomasTram]

Hi Stephane

I have written two functions that should solve your problem, you can find them in the attached .txt file. (I have done some testing, but not extensive, so I apologise in advance for any bugs or weird behaviour.) This is exceptional of course, normally I would just provide some hints to the solution :)

The two functions are
´´´´
int spectra_fast_pk_at_kvec(
struct background * pba,
struct spectra * psp,
double * kvec,
int kvec_size,
double * pk_tot_out /* (must be already allocated with kvec_sizepsp) /
);
''''
and
´´´´
int spectra_fast_pk_at_kvec_and_zvec(
struct background * pba,
struct spectra * psp,
double * kvec,
int kvec_size,
double * zvec,
int zvec_size,
double * pk_tot_out / (must be already allocated with kvec_sizezvec_size) */
''''
The first one returns P(k,z_class) on the grid of log(tau) values found in psp->ln_tau[]. This array is computed by CLASS, and the largest values is slightly before the maximum z required in input. The second function is probably what you want, since it is interpolated in zvec as well. So you just pass the z_array and the k_array you want and it is interpolated.

If you read the code you can see that I made several tricks to increase the interpolation speed. Let me know how it works.

Cheers,
Thomas

PS You need to add the functions to spectra.c and you need to add the functions to spectra.h.

pkfunctions.txt

[plaszczy]

thanks that's great. I will try it and report (this may take some time) stephane Le 31/03/2017 à 10:26, Thomas Tram a écrit :

…

Hi Stephane I have written two functions that should solve your problem, you can find them in the attached .txt file. (I have done some testing, but not extensive, so I apologise in advance for any bugs or weird behaviour.) This is exceptional of course, normally I would just provide some hints to the solution :) The two functions are ´´´´ int spectra_fast_pk_at_kvec( struct background * pba, struct spectra * psp, double * kvec, int kvec_size, double * pk_tot_out /* (must be already allocated with kvec_size/psp) // ); '''' and ´´´´ int spectra_fast_pk_at_kvec_and_zvec( struct background * pba, struct spectra * psp, double * kvec, int kvec_size, double * zvec, int zvec_size, double * pk_tot_out // (must be already allocated with kvec_size/zvec_size) */ '''' The first one returns P(k,z_class) on the grid of log(tau) values found in psp->ln_tau[]. This array is computed by CLASS, and the largest values is slightly before the maximum z required in input. The second function is probably what you want, since it is interpolated in zvec as well. So you just pass the z_array and the k_array you want and it is interpolated. If you read the code you can see that I made several tricks to increase the interpolation speed. Let me know how it works. Cheers, Thomas pkfunctions.txt <https://github.com/lesgourg/class_public/files/884904/pkfunctions.txt> — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#145 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ALVf5iyyHgnPKqjkjdmFob3ftUEeE6KEks5rrLi-gaJpZM4MoFVr>.

[plaszczy]

Hi Thomas, so it compiled smoothly and runs indeed much fatser (about a factor 20 on my bench). But I have troubles accessing the output values that are always 0 (this is maybe just my C++ implementation that use sa std::vector and &vector[0] to access double*, but this generally work) Just to be sure: pk_tot_out should only be sized kvec_size, right? Also a last point: does it worth producing a spectra_fast_pk_at_kvec_at_z with a fixed z values or is it equivalent to using your last function with a z vector of one value? thanks very much stephane Le 31/03/2017 à 10:26, Thomas Tram a écrit :

…

Hi Stephane I have written two functions that should solve your problem, you can find them in the attached .txt file. (I have done some testing, but not extensive, so I apologise in advance for any bugs or weird behaviour.) This is exceptional of course, normally I would just provide some hints to the solution :) The two functions are ´´´´ int spectra_fast_pk_at_kvec( struct background * pba, struct spectra * psp, double * kvec, int kvec_size, double * pk_tot_out /* (must be already allocated with kvec_size/psp) // ); '''' and ´´´´ int spectra_fast_pk_at_kvec_and_zvec( struct background * pba, struct spectra * psp, double * kvec, int kvec_size, double * zvec, int zvec_size, double * pk_tot_out // (must be already allocated with kvec_size/zvec_size) */ '''' The first one returns P(k,z_class) on the grid of log(tau) values found in psp->ln_tau[]. This array is computed by CLASS, and the largest values is slightly before the maximum z required in input. The second function is probably what you want, since it is interpolated in zvec as well. So you just pass the z_array and the k_array you want and it is interpolated. If you read the code you can see that I made several tricks to increase the interpolation speed. Let me know how it works. Cheers, Thomas pkfunctions.txt <https://github.com/lesgourg/class_public/files/884904/pkfunctions.txt> — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#145 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ALVf5iyyHgnPKqjkjdmFob3ftUEeE6KEks5rrLi-gaJpZM4MoFVr>.

[plaszczy]

OK I cured my problem: I was using a z_max_pk =1 in my test (and shoud not with spectra_fast_pk_at_kvec). I guess the output is of size kvecsize*psp->ln_tau_size right?
This one works very well (20 times faster!) and gives the same results than the standard call.
congrats

[ThomasTram]

That's correct. For z_max_pk of 100, psp->ln_tau_size is roughly 450. I think the other function is probably the one you want since you can choose the z-vector yourself. (So this function is potentially much faster if you need less redshifts.)

Cheers,
Thomas

[plaszczy]

OK I confirm the other routine works also well (and is probably less dangerous while as efficient)
I've done a "nl" version (only modifying only the call to spectra_pk_nl_at_z, looks enough could you confirm?

Any chance the code enters the master (or other) branch? it might interest others.
thanks again
stephane

[ThomasTram]

Hi Stephane

Yes that should be fine. I think this new method should be implemented in the master branch in the next release. I will leave this issue open so that I remember!

Cheers,
Thomas

[plaszczy]

Hi Thomas,
I discovered a slight difference between the fast version you provided and the standard one:
below the kmin cutoff (around 1e-5) the former goes directly to 0, while latter is smoothly extrapolated.
This cause problems when one computes integrals with bessel functions (as Cl). So
-is there a way to have some extrapolation in the fast version?
-is there a way to specify (lower) the kmin cutoff?

thanks,
cheers
stephane

[plaszczy]

Hi Thomas,
I don't think actually there is much problem with the kmin cutoff.
It looks I am getting into troubles because I am using CLASS with several threads.
Thread1 computes P(k,z). however thread2 is computing another quantity at z2 that modifies the ba structure so that P(k,z1) gets screwed up. Can you confirm that the function makes indeed use of ba?

I could copy the ba structure at the start of the Pk function so that it does not change during the computation. Is there a way to do this easily?
cheers
stephane