Divide by zero error in project_2d module

[AlexMaraio]

Hi Joe,

I've been using the project_2d module to turn my matter power spectrum P(k) into lensing C_l values and it's generally been working fine, however I've been experiencing some weird errors when I change the distribution of galaxies in redshift bins.
To generate the n(z), I've been using a Gaussian centred at z=1 with some small width, as follows:

redshift_range = np.linspace(0.0, 3.0, 500)
redshift_dist_1 = stats.norm.pdf(redshift_range, loc=1.00, scale=0.15) 

This works fine. However, when I change the Gaussian to be centred at less than z=1, then suddenly the project_2d module starts spitting out this error:

./cosmosis-standard-library/structure/projection/projection_tools/kernel.py:166: RuntimeWarning: divide by zero encountered in divide
  nchi_spline.y * (nchi_spline.x - chi)/nchi_spline.x, 0.)

I'm not really sure what's causing this as I say it's working for higher redshift bins so odd (in my opinion) that it would break when I move to lower redshifts. It also seems to throw off the sampling as I start getting values of Omega_m & sigma_8 from my pipeline that are way off from what they should be (whereas for the z>=1 bins it works fine).

Have you got any ideas what would be causing this strange behaviour?
I can put together a minimum working example if that would help.

Thanks! :)

[joezuntz]

Hi Alex - sorry to just get to this. Can you send me something to reproduce it?

[AlexMaraio]

Hi Joe, sorry it's taken so long for me to get back to this - I've been been doing other things with Cosmosis that didn't get this error, but I'm back to getting this warning again so hence digging out this thread!

I've put together a minimum working example of this error. It's a simple Cosmosis script that generates a Gaussian n(z) with a custom centre, then calls CAMB for the P(k, z), and then pk_to_cl uses both of these (with a trivial likelihood).

In the values.ini file, if I set the z_mean to 1.0 I get no error and it runs fine, however if this is set to a smaller value, such as 0.95, 0.90, or even 0.50, the divide by zero warning appears.

I've also noted that it only appears occasionally, as it seems quite sensitive to the current cosmological parameters that the sampler is processing (hence the seemingly random value of h0).
It seems that the Cl values produced from Pk_to_Cl look okay with no signs of anything obviously breaking (no NaN's or inf's), but still a bit confused by the error being produced.

Hope this example works for you, anything else I can help with just let me know.
DivideByZeroExample.zip

Thanks!

[joezuntz]

Hi Alex,

This doesn't reproduce the warning you mentioned for me. Using the latest version of the cosmosis-standard-library on my laptop it seems to run fine. I've tried randomly varying H0 with the apriori sampler to generate hundreds of different values, but with no change.

What standard library version are you using? It might be something we've fixed (even accidentally!) in project_2d or camb.

Can you try setting debug=T on your machine, to narrow down where this is happening?

Cheers
Joe

[AlexMaraio]

Hi Joe,

I've just downloaded and installed the latest version of the cosmosis-standard-library and still getting the same warning.

It seems like a very odd error, but I think it stems from the chi values in the set_nofchi_splines values function.
It populates the chi values though evaluating the chi_of_z spline over the redshift range, which is a linspace from zero to three in my case.
The minimum value of the chi values seems to change when varying parameters (h0 for my simple case), where I've got values such as 4.40659898e-16 and -8.81319796e-16, and exactly zero. When the minimum value is non-zero, I do not get any divide-by-zero warnings as when we divide by nchi_spline.x, there are no zeros in the array to trigger this warning. However, when the minimum value is zero, the warning is triggered.
A bit concerning that it's giving a negative value for chi at z=0, but it's clearly just a numerical artefact, though I'm surprised this error has not been caught before as I'd have thought that it should always have returned chi=0 at z=0.

The good news is that it doesn't seem to affect the Cl values that it outputs, as it's just slightly changing the very first value for the lensing kernel.
Not sure how best to deal with this. It's quite annoying for this warning to be printed for every sample, so I've put a with np.errstate(divide='ignore', invalid='ignore'): at the start of the get_wofchi_vals function which seems to do the trick. However you might want a more robust solution.

Thanks,
Alex

EDIT: I've been playing around with your Demo 6 of the CFHTLens likelihood, where I've changed it to draw H0 samples at random and managed to get the divide by zero warning to trigger there - so it's not just my code! I'm surprised that this hasn't been raised before?