GSL interpolation error

alanxuhaojie commented 5 years ago

General information

Corrfunc version: latest?
platform: macOX Mojave 10.14.2
installation method (pip/source/other?):source

Issue description

gsl: interp.c:150: ERROR: interpolation error Default GSL error handler invoked. Abort trap: 6 See the following two examples...

Initially, I thought it may be due to my large input redshift (up to ~1). Then I select out the those galaxies with lower z (say, z<0.2) in my samples, it works.

When I am trying to show some examples, it even fails with z<0.2 with random ra, dec. Any idea?

My gsl version is 2.5

Expected behavior

Actual behavior

What have you tried so far?

Minimal failing example

import Corrfunc
[haojiexu] Corrfunc $ ipython
Python 3.7.1 (default, Dec 14 2018, 13:28:58)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.2.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import numpy as np
   ...: import Corrfunc
   ...: from Corrfunc.mocks.DDrppi_mocks import DDrppi_mocks

In [2]: ra=360*np.random.random_sample((50000,))

In [3]: dec=180.*np.random.random_sample((50000,))-90.

In [4]: z = np.random.random_sample((50000,))

In [5]: cz = 300000.*z

In [6]: nthreads=1

In [7]: cosmology=1

In [8]: pimax = 60.0

In [9]: nrpbins = 30
   ...: rpmin = 0.16
   ...: rpmax = 40.0
   ...: bins = np.logspace(np.log10(rpmin), np.log10(rpmax), num=nrpbi
   ...: ns+1, endpoint=True, base=10.0)
   ...: autocorr = 1

In [10]: DD_counts = DDrppi_mocks(autocorr,cosmology,nthreads, pimax,
    ...: bins, ra, dec, cz)
    ...:
gsl: interp.c:150: ERROR: interpolation error
Default GSL error handler invoked.
Abort trap: 6
[haojiexu] Corrfunc $ ipython
Python 3.7.1 (default, Dec 14 2018, 13:28:58)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.2.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import numpy as np
   ...:    ...: import Corrfunc
   ...:    ...: from Corrfunc.mocks.DDrppi_mocks import DDrppi_mocks

In [2]: ra=360*np.random.random_sample((50000,))

In [3]: dec=180.*np.random.random_sample((50000,))-90.

In [4]: z = 0.2*np.random.random_sample((50000,))

In [5]: cz = 300000.*z

In [6]: nthreads=1

In [7]: cosmology=1

In [8]: pimax = 60.0

In [9]: nrpbins = 30
   ...: rpmin = 0.16
   ...: rpmax = 40.0
   ...: bins = np.logspace(np.log10(rpmin), np.log10(rpmax), num=nrpbins+1, endpoint=True, base
   ...: =10.0)
   ...: autocorr = 1
   ...: print('Counting DD pairs...')
   ...: DD_counts = DDrppi_mocks(autocorr,cosmology,nthreads, pimax, bins, ra, dec, cz)
Counting DD pairs...
gsl: interp.c:150: ERROR: interpolation error
Default GSL error handler invoked.
Abort trap: 6
[haojiexu] Corrfunc $

# rest of sample code goes here...

lgarrison commented 5 years ago

Thanks for the report! I can reproduce this bug.

It looks like an error in how we populate the interpolation arrays in set_cosmo_dist(). Specifically, I think this line:

https://github.com/manodeep/Corrfunc/blob/master/utils/set_cosmo_dist.c#L44

should be zmax/max_size instead of 1/max_size. Right now, the interpolation array is effectively capped at 1; we trigger the line 60 break before actually reaching zmax.

Additionally, the interpolation redshifts start at 1/max_size~1e-4, so we will fail the same way if we get a redshift smaller than that. That's causing the error in the zmax=0.2 case.

@manodeep: I think the most robust way to fix this would be to force the last element in the array to be equal to zmax. Probably we should force the first element to be equal to zmin or 0, too

@alanxuhaojie, if you need a quick fix, you can pass comoving distances instead of cz. That will avoid this part of the code entirely.

manodeep commented 5 years ago

@lgarrison Yup - that entire function needs to be updated to use GSL integration. If I remember correctly updating to gsl integration was breaking the tests. That's probably why the gsl integration headers are included in that file

@alanxuhaojie I will second @lgarrison's suggestion. You might be better off simply calculating the co-moving distances with something like astropy coordinates and then set the flag is_comoving_dist=True. Here is the see note in the python wrapper explaining the is_comoving_dist semantics

alanxuhaojie commented 5 years ago

@lgarrison @manodeep I'll go with passing comoving distances instead of cz for now. Thanks for quick responding. Let me know if the bug is fixed.

manodeep / Corrfunc