Correlation function of DR12G SDSS CMASS Catalog

First import all the modules such as healpy and astropy needed for analyzing the structure


In [1]:
import healpix_util as hu
import astropy as ap
import numpy as np
from astropy.io import fits
from astropy.table import Table
import astropy.io.ascii as ascii
from astropy.io import fits
from astropy.constants import c
import matplotlib.pyplot as plt
import math as m
from math import pi
#from scipy.constants import c
import scipy.special as sp
from astroML.decorators import pickle_results
from scipy import integrate
import warnings
from sklearn.neighbors import BallTree
import pickle
import multiprocessing as mp
import time
from cython_metric import *
from lcdmmetric import *
from lcck0metric import *
from progressbar import *
from tqdm import *
from functools import partial
import pymangle
#from astroML.datasets import fetch_sdss_specgals
#from astroML.correlation import bootstrap_two_point_angular
%matplotlib inline

Read the data file (taken from http://cosmo.nyu.edu/~eak306/SDSS-LRG.html ) converted to ascii with comoving distance etc. in V01 reading from pkl files for faster read


In [2]:
data=ascii.read("./output/dr12gcmssrarfLC.dat")

In [3]:
data


Out[3]:
<Table length=230831>
zradecsrardecr
float64float64float64float64float64float64
0.301401321.4803299.973440.2634415.610890.174069
0.458872321.3596919.9333580.3776645.6087850.17337
0.579383322.61418110.1444960.4570345.630680.177055
0.672858325.01105410.5319990.5145345.6725130.183818
0.517716325.20021410.4954980.4172075.6758140.183181
0.485539325.21419210.4176690.3957785.6760580.181823
0.48558325.30530110.5618770.3958055.6776490.18434
0.566184325.42751510.4481610.4486425.6797820.182355
0.51054325.444410.479460.4124675.6800760.182901
0.50439325.46404110.5688650.4083885.6804190.184461
..................
0.521842330.0850592.4922980.4199215.7610710.043499
0.603357317.6646942.6223240.47215.5442950.045768
0.570822317.7920112.7207030.4515995.5465170.047485
0.602232317.9403942.5665820.4713985.5491070.044795
0.570993318.0767132.7061030.4517085.5514860.04723
0.452071318.4410562.5877740.3729915.5578450.045165
0.451214318.4814392.6388980.37245.558550.046057
0.668937318.4161262.7648270.5121875.557410.048255
0.575847318.3835112.7596990.4547935.5568410.048166
0.541101328.6041681.2553760.4324975.7352250.02191

In [4]:
data.remove_column('z')
data.remove_column('ra')
data.remove_column('dec')

In [5]:
data


Out[5]:
<Table length=230831>
srardecr
float64float64float64
0.2634415.610890.174069
0.3776645.6087850.17337
0.4570345.630680.177055
0.5145345.6725130.183818
0.4172075.6758140.183181
0.3957785.6760580.181823
0.3958055.6776490.18434
0.4486425.6797820.182355
0.4124675.6800760.182901
0.4083885.6804190.184461
.........
0.4199215.7610710.043499
0.47215.5442950.045768
0.4515995.5465170.047485
0.4713985.5491070.044795
0.4517085.5514860.04723
0.3729915.5578450.045165
0.37245.558550.046057
0.5121875.557410.048255
0.4547935.5568410.048166
0.4324975.7352250.02191

In [6]:
rs=np.array(data['s'])
rrar=np.array(data['rar'])
rdecr=np.array(data['decr'])

In [7]:
dat=np.array([rs,rrar,rdecr])

In [8]:
dat


Out[8]:
array([[ 0.263441,  0.377664,  0.457034, ...,  0.512187,  0.454793,
         0.432497],
       [ 5.61089 ,  5.608785,  5.63068 , ...,  5.55741 ,  5.556841,
         5.735225],
       [ 0.174069,  0.17337 ,  0.177055, ...,  0.048255,  0.048166,
         0.02191 ]])

In [9]:
dat.reshape(3,len(data))


Out[9]:
array([[ 0.263441,  0.377664,  0.457034, ...,  0.512187,  0.454793,
         0.432497],
       [ 5.61089 ,  5.608785,  5.63068 , ...,  5.55741 ,  5.556841,
         5.735225],
       [ 0.174069,  0.17337 ,  0.177055, ...,  0.048255,  0.048166,
         0.02191 ]])

In [10]:
dat=dat.transpose()

In [11]:
dat


Out[11]:
array([[ 0.263441,  5.61089 ,  0.174069],
       [ 0.377664,  5.608785,  0.17337 ],
       [ 0.457034,  5.63068 ,  0.177055],
       ..., 
       [ 0.512187,  5.55741 ,  0.048255],
       [ 0.454793,  5.556841,  0.048166],
       [ 0.432497,  5.735225,  0.02191 ]])

In [12]:
# Saving the objects:
with open('dr12gcmsLC.pkl', 'w') as f:  # Python 3: open(..., 'wb')
    pickle.dump(dat, f)

In [2]:
# Getting back the objects:
with open('../pkl/dr12gcmsLC.pkl') as f:  # Python 3: open(..., 'rb')
    dat = pickle.load(f)
dat


Out[2]:
array([[ 0.263441,  5.61089 ,  0.174069],
       [ 0.377664,  5.608785,  0.17337 ],
       [ 0.457034,  5.63068 ,  0.177055],
       ..., 
       [ 0.512187,  5.55741 ,  0.048255],
       [ 0.454793,  5.556841,  0.048166],
       [ 0.432497,  5.735225,  0.02191 ]])

In [3]:
LCDMmetricsq(dat[0],dat[1])


Out[3]:
0.013047370022160415

In [4]:
%%time
BT_D = BallTree(dat,metric='pyfunc',func=LCDMmetricsq) 

with open('../pkl/BTDdr12gcmssLCf.pkl', 'w') as f:
    pickle.dump(BT_D,f)


CPU times: user 27.2 s, sys: 260 ms, total: 27.4 s
Wall time: 29.1 s

In [5]:
with open('../pkl/BTDdr12gcmssLCf.pkl') as f:
    BTD = pickle.load(f)
    
BTD


Out[5]:
<sklearn.neighbors.ball_tree.BinaryTree at 0x7f917b622600>

In [6]:
bins=np.arange(0.004,0.084,0.004)
bins


Out[6]:
array([ 0.004,  0.008,  0.012,  0.016,  0.02 ,  0.024,  0.028,  0.032,
        0.036,  0.04 ,  0.044,  0.048,  0.052,  0.056,  0.06 ,  0.064,
        0.068,  0.072,  0.076,  0.08 ])

In [7]:
binsq=bins**2
binsq


Out[7]:
array([  1.60000000e-05,   6.40000000e-05,   1.44000000e-04,
         2.56000000e-04,   4.00000000e-04,   5.76000000e-04,
         7.84000000e-04,   1.02400000e-03,   1.29600000e-03,
         1.60000000e-03,   1.93600000e-03,   2.30400000e-03,
         2.70400000e-03,   3.13600000e-03,   3.60000000e-03,
         4.09600000e-03,   4.62400000e-03,   5.18400000e-03,
         5.77600000e-03,   6.40000000e-03])

In [8]:
%%time
start_time=time.time()
counts_DD=BTD.two_point_correlation(dat,binsq)
print counts_DD
end_time=time.time()
tottime=end_time-start_time
print "Total run time:"
print tottime

with open('../pkl/BTDdr12gcmssDDLCf.pkl', 'w') as f:
    pickle.dump(counts_DD,f)

with open('../pkl/BTDdr12gcmssDDLCf.pkl') as f:
    counts_DD = pickle.load(f)
    
counts_DD


[   1392782    5570266   14834903   31337210   57454827   96059698
  149698703  220491015  310232848  420249246  552033985  706415297
  884578109 1087297679 1315829698 1570443334 1851329932 2159067746
 2492525511 2851367161]
Total run time:
11899.3217258
CPU times: user 3h 15min 50s, sys: 20.6 s, total: 3h 16min 10s
Wall time: 3h 18min 19s

In [9]:
counts_DD


Out[9]:
array([   1392782,    5570266,   14834903,   31337210,   57454827,
         96059698,  149698703,  220491015,  310232848,  420249246,
        552033985,  706415297,  884578109, 1087297679, 1315829698,
       1570443334, 1851329932, 2159067746, 2492525511, 2851367161])

In [10]:
DD=np.diff(counts_DD)

In [11]:
DD


Out[11]:
array([  4177484,   9264637,  16502307,  26117617,  38604871,  53639005,
        70792312,  89741833, 110016398, 131784739, 154381312, 178162812,
       202719570, 228532019, 254613636, 280886598, 307737814, 333457765,
       358841650])

In [12]:
plt.plot(bins[1:len(bins)],DD,'bo-')


Out[12]:
[<matplotlib.lines.Line2D at 0x108a45350>]

In [30]:
RR_zero = (RR == 0)
RR[RR_zero] = 1

In [ ]:


In [ ]:


In [ ]: