Correlation Functions with Nbodykit¶
In [13]:
import numpy as np
import os
import pandas as pd
import seaborn as sns
from astropy.io import fits
# to make this notebook's output stable across runs
np.random.seed(7)
# To plot pretty figures
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
%load_ext autoreload
%autoreload 2
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
In [14]:
DATA_DIR= os.path.join(os.environ['HOME'],
'Downloads/1741p242')
simcat=fits.open(os.path.join(DATA_DIR,
'rs0/obiwan',
'simcat-elg-1741p242.fits'))[1].data
obitractor= fits.open(os.path.join(DATA_DIR,
'rs0/tractor',
'tractor-1741p242.fits'))[1].data
dr5tractor= fits.open(os.path.join(DATA_DIR,
'dr5',
'tractor-1741p242.fits'))[1].data
In [15]:
from astropy.coordinates import SkyCoord
from astropy import units as u
sim_coord = SkyCoord(ra=simcat['ra']*u.degree,
dec=simcat['dec']*u.degree)
obi_coord = SkyCoord(ra=obitractor['ra']*u.degree,
dec=obitractor['dec']*u.degree)
# idx, d2d, d3d = .match_to_catalog_sky(catalog)
# isim,itrac,d= match_radec(simcat.ra, simcat.dec, obitractor.ra, obitractor.dec,
# 1./3600.0,nearest=True)
# not_isim= not_index(isim,len(simcat))
# not_itrac= not_index(itrac,len(obitractor))
In [24]:
sim_coord.match_to_catalog_3d?
Randoms are the Recovered sources¶
In [29]:
isim, itrac, d2d, _ = obi_coord.search_around_sky(sim_coord, u.deg*1./3600)
def indices2bool(indices,n):
keep= np.zeros(n,bool)
keep[indices]=True
return keep
bsim=indices2bool(isim,len(simcat))
btrac=indices2bool(itrac,len(obitractor))
sns.distplot(d2d*3600,color='b')
Out[29]:
<matplotlib.axes._subplots.AxesSubplot at 0x116076c10>
In [30]:
recovered= simcat[bsim]
real= obitractor[~btrac]
Angular Correlation Function¶
In [40]:
## Example
In [45]:
NDATA = 10000
# randomly distributed RA/DEC on the sky
data = RandomCatalog(NDATA, seed=42)
data['RA'] = data.rng.uniform(low=110, high=260, size=len(data))
data['DEC'] = data.rng.uniform(low=-3.6, high=60., size=len(data))
data['Weight'] = np.random.random(size=len(data))
edges = np.linspace(0.1, 10., 20+1) # 20 total bins
# run the algorithm
r_auto = AngularPairCount(data, edges, ra='RA', dec='DEC', weight='Weight')
pc = r_auto.result
plt.plot(pc['theta'], pc['npairs'])
# format the axes
plt.xlabel(r"$\theta$ [$\mathrm{degrees}$]")
plt.ylabel(r"$DD(\theta)$")
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
<ipython-input-45-52ece9314ac0> in <module>()
8
9 # run the algorithm
---> 10 r_auto = AngularPairCount(data, edges, ra='RA', dec='DEC', weight='Weight')
11
12 pc = r_auto.result
/Users/kaylan1/miniconda3/lib/python2.7/site-packages/nbodykit/algorithms/survey_paircount.pyc in __init__(self, first, edges, second, ra, dec, weight, show_progress, **config)
365
366 # run the algorithm
--> 367 self.run()
368
369 def __setstate__(self, state):
/Users/kaylan1/miniconda3/lib/python2.7/site-packages/nbodykit/algorithms/survey_paircount.pyc in run(self)
423 # run
424 sizes = self.comm.allgather(len(pos1))
--> 425 pc = self._run(DDtheta_mocks, kws, sizes, callback=callback)
426
427 # sum results over all ranks
/Users/kaylan1/miniconda3/lib/python2.7/site-packages/nbodykit/algorithms/sim_paircount.pyc in _run(self, func, kwargs, loads, callback)
264 chunks = numpy.array_split(range(loads[self.comm.rank]), N, axis=0)
265 for i, chunk in enumerate(chunks):
--> 266 this_pc = run(chunk)
267 if self.comm.rank == largest_load:
268 self.logger.info("%d%% done" % (N*(i+1)))
/Users/kaylan1/miniconda3/lib/python2.7/site-packages/nbodykit/algorithms/sim_paircount.pyc in run(chunk)
250 if callback is not None:
251 callback(kwargs, chunk)
--> 252 return self._run_corrfunc(func, **kwargs)
253
254 # log the function start
/Users/kaylan1/miniconda3/lib/python2.7/site-packages/nbodykit/algorithms/sim_paircount.pyc in _run_corrfunc(self, func, **kws)
298 msg += "stdout: %s\n" % stdout
299 msg += "stderr: %s" % stderr
--> 300 raise RuntimeError(msg)
301
302 return CorrfuncResult(result)
RuntimeError: calling the function 'Corrfunc.mocks.DDtheta_mocks.DDtheta_mocks' failed:
exception: Could not import the C extension for the angular correlation function for mocks.
stdout:
stderr:
In [41]:
from nbodykit.source.catalog.array import ArrayCatalog,RandomCatalog
recovered= ArrayCatalog(simcat[bsim])
real= ArrayCatalog(obitractor[~btrac])
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
<ipython-input-41-acde3c3842fc> in <module>()
----> 1 from nbodykit.source.catalog.array import ArrayCatalog,RandomCatalog
2 recovered= ArrayCatalog(simcat[bsim])
3 real= ArrayCatalog(obitractor[~btrac])
ImportError: cannot import name RandomCatalog
In [42]:
# compute the cross correlation pair counts
from nbodykit.lab import AngularPairCount,RandomCatalog
edges = np.linspace(0.25/1e2,0.25/1e1,20+1)
r_cross = AngularPairCount(real, edges, ra='ra', dec='dec',
second=recovered)
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
<ipython-input-42-a66451d4a86e> in <module>()
3 edges = np.linspace(0.25/1e2,0.25/1e1,20+1)
4 r_cross = AngularPairCount(real, edges, ra='ra', dec='dec',
----> 5 second=recovered)
/Users/kaylan1/miniconda3/lib/python2.7/site-packages/nbodykit/algorithms/survey_paircount.pyc in __init__(self, first, edges, second, ra, dec, weight, show_progress, **config)
365
366 # run the algorithm
--> 367 self.run()
368
369 def __setstate__(self, state):
/Users/kaylan1/miniconda3/lib/python2.7/site-packages/nbodykit/algorithms/survey_paircount.pyc in run(self)
423 # run
424 sizes = self.comm.allgather(len(pos1))
--> 425 pc = self._run(DDtheta_mocks, kws, sizes, callback=callback)
426
427 # sum results over all ranks
/Users/kaylan1/miniconda3/lib/python2.7/site-packages/nbodykit/algorithms/sim_paircount.pyc in _run(self, func, kwargs, loads, callback)
264 chunks = numpy.array_split(range(loads[self.comm.rank]), N, axis=0)
265 for i, chunk in enumerate(chunks):
--> 266 this_pc = run(chunk)
267 if self.comm.rank == largest_load:
268 self.logger.info("%d%% done" % (N*(i+1)))
/Users/kaylan1/miniconda3/lib/python2.7/site-packages/nbodykit/algorithms/sim_paircount.pyc in run(chunk)
250 if callback is not None:
251 callback(kwargs, chunk)
--> 252 return self._run_corrfunc(func, **kwargs)
253
254 # log the function start
/Users/kaylan1/miniconda3/lib/python2.7/site-packages/nbodykit/algorithms/sim_paircount.pyc in _run_corrfunc(self, func, **kws)
298 msg += "stdout: %s\n" % stdout
299 msg += "stderr: %s" % stderr
--> 300 raise RuntimeError(msg)
301
302 return CorrfuncResult(result)
RuntimeError: calling the function 'Corrfunc.mocks.DDtheta_mocks.DDtheta_mocks' failed:
exception: Could not import the C extension for the angular correlation function for mocks.
stdout:
stderr: