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
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 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]:
# Getting back the objects:
with open('../output/datsLCDM.pkl') as f: # Python 3: open(..., 'rb')
dat = pickle.load(f)
dat
Out[2]:
In [3]:
bins=np.arange(0.001,0.0801,0.001)
In [4]:
print bins
In [5]:
binsq=bins**2
binsq
Out[5]:
In [9]:
%%time
BT_D = BallTree(dat,metric='pyfunc',func=LCDMmetricsq,leaf_size=5)
with open('../pkl/BTDdatsLCDMsq.pkl', 'w') as f:
pickle.dump(BT_D,f)
with open('../pkl/BTDdatsLCDMsq.pkl') as f:
BTD = pickle.load(f)
BTD
In [11]:
%%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/BTDR72sDDLCDM.pkl', 'w') as f:
pickle.dump(counts_DD,f)
In [12]:
with open('../pkl/BTDR72sDDLCDM.pkl') as f:
counts_DD = pickle.load(f)
counts_DD
Out[12]:
In [14]:
DD=np.diff(counts_DD)
DD
Out[14]:
In [15]:
len(dat)
Out[15]:
BallTree.two_point_correlation works almost 10 times faster! with leaf_size=5 Going with it to the random catalog
In [16]:
# Getting back the objects:
with open('../output/datRs100kLCDM.pkl') as f: # Python 3: open(..., 'rb')
datR = pickle.load(f)
datR
Out[16]:
In [17]:
%%time
BT_R = BallTree(dat,metric='pyfunc',func=LCDMmetricsq,leaf_size=5)
with open('../pkl/BTRdatsLCDMsq.pkl', 'w') as f:
pickle.dump(BT_R,f)
with open('../pkl/BTRdatsLCDMsq.pkl') as f:
BTR = pickle.load(f)
BTR
In [18]:
%%time
start_time=time.time()
counts_RR=BTR.two_point_correlation(datR,binsq)
print counts_RR
end_time=time.time()
tottime=end_time-start_time
print "Total run time:"
print tottime
with open('../pkl/BTRDR72100kRRLCDMsq.pkl', 'w') as f:
pickle.dump(counts_RR,f)
with open('../pkl/BTRDR72100kRRLCDMsq.pkl') as f:
counts_RR = pickle.load(f)
counts_RR
In [19]:
RR=np.diff(counts_RR)
RR
Out[19]:
In [41]:
%%time
start_time=time.time()
counts_DR=BTD.two_point_correlation(datR,binsq)
end_time=time.time()
tottime=end_time-start_time
print "Total run time:"
print tottime
with open('../pkl/BTRDR72100kDRLCDMsq_.pkl', 'w') as f:
pickle.dump(counts_DR,f)
In [42]:
with open('../pkl/BTRDR72100kDRLCDMsq_.pkl') as f:
counts_DR = pickle.load(f)
counts_DR
Out[42]:
In [43]:
DR=np.diff(counts_DR)
DR
Out[43]:
In [44]:
correl100k=(DD-2.0*DR+RR)/RR
correl100k
Out[44]:
In [45]:
plt.plot(bins[1:len(bins)],correl100k)
Out[45]:
In [47]:
plt.yscale('log')
plt.plot(c*1e-5*bins[1:len(bins)],correl100k)
Out[47]:
In [48]:
plt.xscale('log')
plt.yscale('log')
plt.plot(c*1e-5*bins[1:len(bins)],correl100k)
Out[48]:
In [32]:
(DD*1.0/RR)
Out[32]:
In [35]:
(DR*1.0)/RR
Out[35]:
In [39]:
plt.plot(bins[1:len(bins)],(DR*1.0)/RR-1.0)
Out[39]:
In [40]:
#plt.xscale('log')
plt.yscale('log')
plt.plot(c*1e-5*bins[1:len(bins)],(DR*1.0)/RR-1.0)
Out[40]:
In [ ]:
In [36]:
DD-DR
Out[36]:
In [37]:
DD
Out[37]:
In [38]:
DR
Out[38]:
In [ ]:
In [49]:
dataR=ascii.read("/Users/rohin/Downloads/random-DR7-Full.ascii")
In [50]:
Ez = lambda x: 1/m.sqrt(0.3*(1+x)**3+0.7)
np.vectorize(Ez)
#Calculate comoving distance of a data point using the Redshift - This definition is based on the cosmology model we take. Here the distance for E-dS universe is considered. Also note that c/H0 ratio is cancelled in the equations and hence not taken.
def DC_LCDM(z):
return integrate.quad(Ez, 0, z)[0]
DC_LCDM=np.vectorize(DC_LCDM)
dr12gcmn = open("../output/rdr12gcmnsrarf.dat",'w')
dr12gcmn.write("z\t ra\t dec\t s\t rar\t decr \n")
for i in range(0,len(zr)):
dr12gcmn.write("%f\t " %zr[i])
dr12gcmn.write("%f\t %f\t " %(rar[i],decr[i]))
dr12gcmn.write("%f\t " %DC_LCDM(zr[i]))
dr12gcmn.write("%f\t %f\n " %(rar[i]*pi/180.0,decr[i]*pi/180.0))
dr12gcmn.close()
Out[50]:
In [51]:
import pyfits
dataRf=pyfits.open("/Users/rohin/Downloads/random-DR7-Full.fits")
dataRf=dataRf[1].data
dataRf.columns
Out[51]:
In [52]:
z=dataRf['Z']
ra=dataRf['RA']
dec=dataRf['DEC']
In [53]:
Ez = lambda x: 1/m.sqrt(0.3*(1+x)**3+0.7)
np.vectorize(Ez)
#Calculate comoving distance of a data point using the Redshift - This definition is based on the cosmology model we take. Here the distance for E-dS universe is considered. Also note that c/H0 ratio is cancelled in the equations and hence not taken.
def DC_LCDM(z):
return integrate.quad(Ez, 0, z)[0]
DC_LCDM=np.vectorize(DC_LCDM)
In [54]:
dr72r = open("../output/DR72srandf16M.dat",'w')
dr72r.write("z\t ra\t dec\t s\t rar\t decr \n")
for i in range(0,len(z)):
dr72r.write("%f\t " %z[i])
dr72r.write("%f\t %f\t " %(ra[i],dec[i]))
dr72r.write("%f\t " %DC_LCDM(z[i]))
dr72r.write("%f\t %f\n " %(ra[i]*pi/180.0,dec[i]*pi/180.0))
dr72r.close()
In [55]:
z
Out[55]:
In [6]:
dataR=ascii.read("../output/DR72srandf16M.dat")
In [7]:
dataR
Out[7]:
In [8]:
dataR.remove_column('z')
dataR.remove_column('ra')
dataR.remove_column('dec')
In [9]:
dataR
Out[9]:
In [14]:
rs=np.array(dataR['s'])
rrar=np.array(dataR['rar'])
rdecr=np.array(dataR['decr'])
In [15]:
datR=np.array([rs,rrar,rdecr])
In [16]:
datR
Out[16]:
In [17]:
datR.reshape(3,len(dataR))
Out[17]:
In [18]:
datR=datR.transpose()
In [19]:
datR
Out[19]:
In [66]:
# Saving the objects:
with open('../pkl/datRs16MLCDM.pkl', 'w') as f: # Python 3: open(..., 'wb')
pickle.dump(datR, f)
In [6]:
with open('../pkl/datRs16MLCDM.pkl') as f: # Python 3: open(..., 'wb')
datR=pickle.load(f)
In [7]:
len(datR)
Out[7]:
In [20]:
%%time
BTR = BallTree(datR,metric='pyfunc',func=LCDMmetricsq,leaf_size=5)
In [ ]:
with open('../pkl/BTR16MdatsLCDM.pkl', 'w') as f:
pickle.dump(BT_R2,f)
In [8]:
with open('../pkl/BTR16MdatsLCDM.pkl') as f:
BTR = pickle.load(f)
BTR
Out[8]:
In [ ]:
%%time
start_time=time.time()
counts_RR=BTR.two_point_correlation(datR,binsq)
print counts_RR
end_time=time.time()
tottime=end_time-start_time
print "Total run time:"
print tottime
with open('../pkl/BTR16MRRLCDM.pkl', 'w') as f:
pickle.dump(counts_RR,f)
In [22]:
with open('../pkl/BTR16MRRLCDM.pkl') as f:
counts_RR = pickle.load(f)
counts_RR
Out[22]:
In [23]:
RR=np.diff(counts_RR)
RR
Out[23]:
In [ ]:
%%time
start_time=time.time()
counts_DR=BTR.two_point_correlation(dat,binsq)
print counts_DR
end_time=time.time()
tottime=end_time-start_time
print "Total run time:"
print tottime
with open('../pkl/BTR16MDRLCDM.pkl', 'w') as f:
pickle.dump(counts_DR,f)
In [24]:
with open('../pkl/BTR16MDRLCDM.pkl') as f:
counts_DR = pickle.load(f)
counts_DR
Out[24]:
In [25]:
DR=np.diff(counts_DR)
DR
Out[25]:
In [26]:
with open('../pkl/BTDR72sDDLCDM.pkl') as f:
counts_DD = pickle.load(f)
counts_DD
Out[26]:
In [27]:
DD=np.diff(counts_DD)
DD
Out[27]:
In [37]:
len(dataR)
Out[37]:
In [38]:
len(dat)
Out[38]:
In [39]:
NrbyNd=1.0*len(dataR)/len(dat)
NrbyNd
Out[39]:
In [30]:
corrells16m=(DD*1.0-2.0*DR+RR*1.0)/RR
In [31]:
corrells16m
Out[31]:
In [33]:
plt.plot(bins[1:len(bins)],corrells16m,'bo-')
Out[33]:
In [35]:
plt.yscale('log')
plt.plot(bins[1:len(bins)],corrells16m,'b.-')
Out[35]:
In [36]:
plt.xscale('log')
plt.yscale('log')
plt.plot(c*1e-5*bins[1:len(bins)],corrells16m,'b.-')
Out[36]:
In [42]:
corrells16mr=(DD*NrbyNd**2-2.0*NrbyNd*DR+RR*1.0)/RR
In [43]:
corrells16mr
Out[43]:
In [ ]:
In [46]:
plt.yscale('log')
plt.plot(bins[1:len(bins)],(1.0+corrells16mr),'b.-')
Out[46]:
In [ ]:
In [ ]:
In [ ]:
In [53]:
RR_zero = (RR == 0)
RR[RR_zero] = 1
In [54]:
correl=1.0*DD/RR-1.0
In [55]:
correl
Out[55]:
In [64]:
plt.plot(bins[1:len(bins)],correl,'go-')
Out[64]:
In [67]:
plt.yscale('log')
plt.plot(bins[1:len(bins)],correl,'go-')
Out[67]:
In [63]:
plt.plot(bins[1:len(bins)],bins[1:len(bins)]**2*correl,'go-')
Out[63]:
In [58]:
plt.plot(bins[2:len(bins)],bins[2:len(bins)]*bins[2:len(bins)]*correl[1:len(bins)],'go')
Out[58]:
In [59]:
plt.plot(bins[2:len(bins)],correl[1:len(bins)],'go')
Out[59]:
In [60]:
plt.plot(bins[2:len(bins)],correl[1:len(bins)],'go-')
plt.savefig("correl1.pdf")
In [61]:
plt.plot(bins[2:len(bins)]*299792.458/100,correl[1:len(bins)],'bo-')
Out[61]:
In [ ]:
In [ ]:
In [ ]:
start_time=time.time()
arg=[(dat,bins)]
pool=mp.Pool(8)
#%timeit
#@pickle_results("DR72DD2ptc.pkl")
def mf_wrap(args):
return BTD.two_point_correlation(*args)
%timeit counts_DD=pool.map(mf_wrap,arg)
end_time=time.time()
tottime=end_time-start_time
print "Total run time:"
print tottime
In [ ]:
from functools import partial
def harvester(text, case):
X = case[0]
return text + str(X)
partial_harvester = partial(harvester, case=RAW_DATASET)
partial_qr=partial(BTD.query_radius,count_only=True)
if __name__ == '__main__':
pool = multiprocessing.Pool(processes=6)
case_data = RAW_DATASET
pool.map(partial_harvester, case_data, 1)
pool.close()
pool.join()
mapfunc = partial(BTD.query_radius, count_only=True)
map(mapfunc, volume_ids)
In [ ]:
#ascii.write("DR72DDbinned.dat",(bins[1:len(bins)],DDresult))
start_time=time.time()
@pickle_results("DR72DDmp1.pkl")
def ddcal(BTD,dat,bins,Nbins):
counts_DD=np.zeros(Nbins)
for i in tqdm(range(Nbins)):
counts_DD[i]=np.sum(BTD.query_radius(dat, bins[i],count_only=True))
DD = np.diff(counts_DD)
print counts_DD
print DD
return DD
def mf_wrap(args):
return ddcal(*args)
pool=mp.Pool(8)
arg=[(BTD,dat,bins,Nbins)]
%timeit DDresult=pool.map(mf_wrap,arg)
#DDresult = ddcal(BTD,dat,bins,Nbins)
end_time=time.time()
tottime=end_time-start_time
print "Total run time:"
print tottime
In [ ]:
%timeit dat
In [ ]:
DDresult[0]
In [ ]:
DDresult[1]
In [ ]:
plt.plot(bins[1:len(bins)],DDresult[0],'ro')
In [ ]:
In [ ]:
In [ ]:
def myfun(a,b):
print a + b
return a+b
def mf_wrap(args):
return myfun(*args)
p = mp.Pool(4)
fl = [(a,b) for a in range(3) for b in range(2)]
p.map(mf_wrap, fl)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
counts_DD=np.zeros(Nbins)
for i in range(Nbins):
counts_DD[i]=np.sum(BTD.query_radius(dat, bins[i],count_only=True))
DD = np.diff(counts_DD)
In [ ]:
print counts_DD
print DD
In [ ]:
plt.plot(bins[1:len(bins)],DD,'ro')
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
dataR=fits.open("/Users/rohin/Downloads/random-DR7-Full.fits")
In [ ]:
dataR=dataR[1].data
In [ ]:
len(dataR)
In [ ]:
In [ ]:
tdata=np.array(data)
In [ ]:
type(tdata[4])
In [ ]:
tdata.shape
In [ ]:
In [ ]:
tdata.shape
In [ ]:
tdata=np.atleast_d(tdata)
In [ ]:
tdata.shape
In [ ]:
tdata.reshape(len(tdata),3)
In [ ]:
tdata=np.asarray(data)
tdata=tdata.transpose()
In [ ]:
In [ ]:
In [ ]:
tdata
In [ ]:
len(tdata)
In [ ]:
In [ ]:
BTD.two_point_correlationpoint_correlationpoint_correlationpoint_correlationtime
stime=time.time()
tpcf=BTD.two_point_correlation(dat,bins)
print time.time()-stime
print tpcf
plt.plot(bins,tpcf)
In [ ]:
stime=time.time()
tpcfd=BTD.two_point_correlation(dat,bins,dualtree=True)
print time.time()-stime
print tpcfd
plt.plot(bins,tpcfd)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
X
In [ ]:
In [ ]:
np.random.seed(0)
X = np.random.random((30,3))
r = np.linspace(0, 1, 10)
tree = BallTree(X,metric='pyfunc',func=LCDMmetric)
s = pickle.dumps(tree)
treedump = pickle.loads(s)
treedump.two_point_correlation(X,r)
In [ ]:
BT_D = BallTree(data)
BT_R = BallTree(data_R)
counts_DD = np.zeros(Nbins + 1)
counts_RR = np.zeros(Nbins + 1)
for i in range(Nbins + 1):
counts_DD[i] = np.sum(BT_D.query_radius(data, bins[i],
count_only=True))
counts_RR[i] = np.sum(BT_R.query_radius(data_R, bins[i],
count_only=True))
DD = np.diff(counts_DD)
RR = np.diff(counts_RR)
# check for zero in the denominator
RR_zero = (RR == 0)
RR[RR_zero] = 1
if method == 'standard':
corr = factor ** 2 * DD / RR - 1
elif method == 'landy-szalay':
if sklearn_has_two_point:
counts_DR = KDT_R.two_point_correlation(data, bins)
else:
counts_DR = np.zeros(Nbins + 1)
for i in range(Nbins + 1):
counts_DR[i] = np.sum(BT_R.query_radius(data, bins[i],
count_only=True))
DR = np.diff(counts_DR)
corr = (factor ** 2 * DD - 2 * factor * DR + RR) / RR
corr[RR_zero] = np.nan
return corr
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
dr7fdat=np.array([data['s'][0:300] data['rar'][0:300] data['decr'][0:300]])
dr7fdat
In [ ]:
dr7fdat[2]
In [ ]:
In [ ]:
def LCDMmetric(p1,p2):
costheta=m.sin(dec1)*m.sin(dec2)+m.cos(dec1)*m.cos(dec2)*m.cos(ra1-ra2)
s1=DC_LCDM(z1)
s2=DC_LCDM(z2)
return np.sqrt(s1**2+s2**2-2.0*s1*s2*costheta)
In [ ]:
In [ ]:
#fdata=fits.open("/Users/rohin/Downloads/DR7-Full.fits")
In [ ]:
#fdata.writeto("./output/DR7fulltrim.fits")
In [ ]:
fdata=fits.open("./output/DR7fulltrim.fits")
In [ ]:
cols=fdata[1].columns
In [ ]:
cols.del_col('ZTYPE')
In [ ]:
cols.del_col('SECTOR')
cols.del_col('FGOTMAIN')
cols.del_col('QUALITY')
cols.del_col('ISBAD')
cols.del_col('M')
cols.del_col('MMAX')
cols.del_col('ILSS')
cols.del_col('ICOMB')
cols.del_col('VAGC_SELECT')
cols.del_col('LSS_INDEX')
cols.del_col('FIBERWEIGHT')
cols.del_col('PRIMTARGET')
cols.del_col('MG')
cols.del_col('SECTOR_COMPLETENESS')
cols.del_col('COMOV_DENSITY')
cols.del_col('RADIAL_WEIGHT')
In [ ]:
fdata[1].columns
In [ ]:
fdata.writeto("./output/DR7fullzradec.fits")
In [ ]:
fdat=fits.open("./output/DR7fullzradec.fits")
In [ ]:
fdat[1].columns
In [ ]:
fdat[1].data['Z']
In [ ]:
fdat[1].data['RA']
In [ ]:
In [ ]:
comovlcdm=DC_LCDM(fdat[1].data['Z'])
In [ ]:
fdat[1].data['Z']
In [ ]:
comovlcdm
In [ ]:
comovlcdm.dtype
In [ ]:
#cols=fdat[1].columns
In [ ]:
nc=fits.Column(name='COMOV',format='D',array=comovlcdm)
In [ ]:
nc1=fits.Column(name='COMOV',format='D')
In [ ]:
fdata[1].data['Z']
In [ ]:
fdata[1].data['RA']
In [ ]:
nc
In [ ]:
nc.dtype
In [ ]:
#cols.add_col(nc)
In [ ]:
fdat[1].columns
In [ ]:
fdat[1].columns.info()
In [ ]:
fdat[1].columns.add_col(nc1)
In [ ]:
fdat[1].data['COMOV']=comovlcdm
In [ ]:
comovlcdm
In [ ]:
fdat[1].data['Z']
In [ ]:
fdat[1].data['COMOV']
In [ ]:
fdat[1].data['RA']
In [ ]:
fdat[1].data['RA']=fdat[1].data['RA']*pi/180.0
In [ ]:
comovlcdm=DC_LCDM(fdat[1].data['Z'])
comovlcdm
In [ ]:
In [ ]:
Random catalog created based on the survey limitations also taken from http://cosmo.nyu.edu/~eak306/SDSS-LRG.html
In [ ]:
dataR=fits.open("/Users/rohin/Downloads/random-DR7-Full.fits")
In [ ]:
dataR
In [ ]:
dataR=dataR[1].data
In [ ]:
len(dataR)
In [ ]:
NSIDE=512
dr72hpix=hu.HealPix("ring",NSIDE)
In [ ]:
pixdata = open("./output/pixdatadr72VAGCfullrand.dat",'w')
pixdata.write("z\t pix \n")
for i in range(0,len(data)-1):
pixdata.write("%f\t" %data['z'][i])
pixdata.write("%d\n" %dr72hpix.eq2pix(dataR['ra'][i],dataR['dec'][i]))
pixdata.close()
In [ ]:
pixdata = ascii.read("./output/pixdatadr72VAGCfullrand.dat")
hpixdata=np.array(np.zeros(hu.nside2npix(NSIDE)))
for j in range(len(pixdata)):
hpixdata[pixdata[j]['pix']]+=1
In [ ]:
hpixdata
In [ ]:
hu.mollview(hpixdata,rot=180)
In [ ]:
hu.orthview(hpixdata)
In [ ]:
"""
Tools for computing two-point correlation functions.
"""
#from .utils import check_random_state
# From scikit-learn utilities:
def check_random_state(seed):
"""Turn seed into a np.random.RandomState instance
If seed is None, return the RandomState singleton used by np.random.
If seed is an int, return a new RandomState instance seeded with seed.
If seed is already a RandomState instance, return it.
Otherwise raise ValueError.
"""
if seed is None or seed is np.random:
return np.random.mtrand._rand
if isinstance(seed, (int, np.integer)):
return np.random.RandomState(seed)
if isinstance(seed, np.random.RandomState):
return seed
raise ValueError('%r cannot be used to seed a numpy.random.RandomState'
' instance' % seed)
# Check if scikit-learn's two-point functionality is available.
# This was added in scikit-learn version 0.14
try:
from sklearn.neighbors import KDTree
sklearn_has_two_point = True
except ImportError:
import warnings
sklearn_has_two_point = False
def uniform_sphere(RAlim, DEClim, size=1):
"""Draw a uniform sample on a sphere
Parameters
----------
RAlim : tuple
select Right Ascension between RAlim[0] and RAlim[1]
units are degrees
DEClim : tuple
select Declination between DEClim[0] and DEClim[1]
size : int (optional)
the size of the random arrays to return (default = 1)
Returns
-------
RA, DEC : ndarray
the random sample on the sphere within the given limits.
arrays have shape equal to size.
"""
zlim = np.sin(np.pi * np.asarray(DEClim) / 180.)
z = zlim[0] + (zlim[1] - zlim[0]) * np.random.random(size)
DEC = (180. / np.pi) * np.arcsin(z)
RA = RAlim[0] + (RAlim[1] - RAlim[0]) * np.random.random(size)
return RA, DEC
def ra_dec_to_xyz(ra, dec):
"""Convert ra & dec to Euclidean points
Parameters
----------
ra, dec : ndarrays
Returns
x, y, z : ndarrays
"""
sin_ra = np.sin(ra * np.pi / 180.)
cos_ra = np.cos(ra * np.pi / 180.)
sin_dec = np.sin(np.pi / 2. - dec * np.pi / 180.)
cos_dec = np.cos(np.pi / 2. - dec * np.pi / 180.)
return (cos_ra * sin_dec,
sin_ra * sin_dec,
cos_dec)
def angular_dist_to_euclidean_dist(D, r=1):
"""convert angular distances to euclidean distances"""
return 2 * r * np.sin(0.5 * D * np.pi / 180.)
def two_point(data, bins, method='standard',
data_R=None, random_state=None):
"""Two-point correlation function
Parameters
----------
data : array_like
input data, shape = [n_samples, n_features]
bins : array_like
bins within which to compute the 2-point correlation.
shape = Nbins + 1
method : string
"standard" or "landy-szalay".
data_R : array_like (optional)
if specified, use this as the random comparison sample
random_state : integer, np.random.RandomState, or None
specify the random state to use for generating background
Returns
-------
corr : ndarray
the estimate of the correlation function within each bin
shape = Nbins
"""
data = np.asarray(data)
bins = np.asarray(bins)
rng = check_random_state(random_state)
if method not in ['standard', 'landy-szalay']:
raise ValueError("method must be 'standard' or 'landy-szalay'")
if bins.ndim != 1:
raise ValueError("bins must be a 1D array")
if data.ndim == 1:
data = data[:, np.newaxis]
elif data.ndim != 2:
raise ValueError("data should be 1D or 2D")
n_samples, n_features = data.shape
Nbins = len(bins) - 1
# shuffle all but one axis to get background distribution
if data_R is None:
data_R = data.copy()
for i in range(n_features - 1):
rng.shuffle(data_R[:, i])
else:
data_R = np.asarray(data_R)
if (data_R.ndim != 2) or (data_R.shape[-1] != n_features):
raise ValueError('data_R must have same n_features as data')
factor = len(data_R) * 1. / len(data)
if sklearn_has_two_point:
# Fast two-point correlation functions added in scikit-learn v. 0.14
KDT_D = KDTree(data)
KDT_R = KDTree(data_R)
counts_DD = KDT_D.two_point_correlation(data, bins)
counts_RR = KDT_R.two_point_correlation(data_R, bins)
else:
warnings.warn("Version 0.3 of astroML will require scikit-learn "
"version 0.14 or higher for correlation function "
"calculations. Upgrade to sklearn 0.14+ now for much "
"faster correlation function calculations.")
BT_D = BallTree(data)
BT_R = BallTree(data_R)
counts_DD = np.zeros(Nbins + 1)
counts_RR = np.zeros(Nbins + 1)
for i in range(Nbins + 1):
counts_DD[i] = np.sum(BT_D.query_radius(data, bins[i],
count_only=True))
counts_RR[i] = np.sum(BT_R.query_radius(data_R, bins[i],
count_only=True))
DD = np.diff(counts_DD)
RR = np.diff(counts_RR)
# check for zero in the denominator
RR_zero = (RR == 0)
RR[RR_zero] = 1
if method == 'standard':
corr = factor ** 2 * DD / RR - 1
elif method == 'landy-szalay':
if sklearn_has_two_point:
counts_DR = KDT_R.two_point_correlation(data, bins)
else:
counts_DR = np.zeros(Nbins + 1)
for i in range(Nbins + 1):
counts_DR[i] = np.sum(BT_R.query_radius(data, bins[i],
count_only=True))
DR = np.diff(counts_DR)
corr = (factor ** 2 * DD - 2 * factor * DR + RR) / RR
corr[RR_zero] = np.nan
return corr
def bootstrap_two_point(data, bins, Nbootstrap=10,
method='standard', return_bootstraps=False,
random_state=None):
"""Bootstrapped two-point correlation function
Parameters
----------
data : array_like
input data, shape = [n_samples, n_features]
bins : array_like
bins within which to compute the 2-point correlation.
shape = Nbins + 1
Nbootstrap : integer
number of bootstrap resamples to perform (default = 10)
method : string
"standard" or "landy-szalay".
return_bootstraps: bool
if True, return full bootstrapped samples
random_state : integer, np.random.RandomState, or None
specify the random state to use for generating background
Returns
-------
corr, corr_err : ndarrays
the estimate of the correlation function and the bootstrap
error within each bin. shape = Nbins
"""
data = np.asarray(data)
bins = np.asarray(bins)
rng = check_random_state(random_state)
if method not in ['standard', 'landy-szalay']:
raise ValueError("method must be 'standard' or 'landy-szalay'")
if bins.ndim != 1:
raise ValueError("bins must be a 1D array")
if data.ndim == 1:
data = data[:, np.newaxis]
elif data.ndim != 2:
raise ValueError("data should be 1D or 2D")
if Nbootstrap < 2:
raise ValueError("Nbootstrap must be greater than 1")
n_samples, n_features = data.shape
# get the baseline estimate
corr = two_point(data, bins, method=method, random_state=rng)
bootstraps = np.zeros((Nbootstrap, len(corr)))
for i in range(Nbootstrap):
indices = rng.randint(0, n_samples, n_samples)
bootstraps[i] = two_point(data[indices, :], bins, method=method,
random_state=rng)
# use masked std dev in case of NaNs
corr_err = np.asarray(np.ma.masked_invalid(bootstraps).std(0, ddof=1))
if return_bootstraps:
return corr, corr_err, bootstraps
else:
return corr, corr_err
def two_point_angular(ra, dec, bins, method='standard', random_state=None):
"""Angular two-point correlation function
A separate function is needed because angular distances are not
euclidean, and random sampling needs to take into account the
spherical volume element.
Parameters
----------
ra : array_like
input right ascention, shape = (n_samples,)
dec : array_like
input declination
bins : array_like
bins within which to compute the 2-point correlation.
shape = Nbins + 1
method : string
"standard" or "landy-szalay".
random_state : integer, np.random.RandomState, or None
specify the random state to use for generating background
Returns
-------
corr : ndarray
the estimate of the correlation function within each bin
shape = Nbins
"""
ra = np.asarray(ra)
dec = np.asarray(dec)
rng = check_random_state(random_state)
if method not in ['standard', 'landy-szalay']:
raise ValueError("method must be 'standard' or 'landy-szalay'")
if bins.ndim != 1:
raise ValueError("bins must be a 1D array")
if (ra.ndim != 1) or (dec.ndim != 1) or (ra.shape != dec.shape):
raise ValueError('ra and dec must be 1-dimensional '
'arrays of the same length')
n_features = len(ra)
Nbins = len(bins) - 1
# draw a random sample with N points
ra_R, dec_R = uniform_sphere((min(ra), max(ra)),
(min(dec), max(dec)),
2 * len(ra))
data = np.asarray(ra_dec_to_xyz(ra, dec), order='F').T
data_R = np.asarray(ra_dec_to_xyz(ra_R, dec_R), order='F').T
# convert spherical bins to cartesian bins
bins_transform = angular_dist_to_euclidean_dist(bins)
return two_point(data, bins_transform, method=method,
data_R=data_R, random_state=rng)
def bootstrap_two_point_angular(ra, dec, bins, method='standard',
Nbootstraps=10, random_state=None):
# type: (object, object, object, object, object, object) -> object
"""Angular two-point correlation function
A separate function is needed because angular distances are not
euclidean, and random sampling needs to take into account the
spherical volume element.
Parameters
----------
ra : array_like
input right ascention, shape = (n_samples,)
dec : array_like
input declination
bins : array_like
bins within which to compute the 2-point correlation.
shape = Nbins + 1
method : string
"standard" or "landy-szalay".
Nbootstraps : int
number of bootstrap resamples
random_state : integer, np.random.RandomState, or None
specify the random state to use for generating background
Returns
-------
corr : ndarray
the estimate of the correlation function within each bin
shape = Nbins
dcorr : ndarray
error estimate on dcorr (sample standard deviation of
bootstrap resamples)
bootstraps : ndarray
The full sample of bootstraps used to compute corr and dcorr
"""
ra = np.asarray(ra)
dec = np.asarray(dec)
rng = check_random_state(random_state)
if method not in ['standard', 'landy-szalay']:
raise ValueError("method must be 'standard' or 'landy-szalay'")
if bins.ndim != 1:
raise ValueError("bins must be a 1D array")
if (ra.ndim != 1) or (dec.ndim != 1) or (ra.shape != dec.shape):
raise ValueError('ra and dec must be 1-dimensional '
'arrays of the same length')
n_features = len(ra)
Nbins = len(bins) - 1
data = np.asarray(ra_dec_to_xyz(ra, dec), order='F').T
# convert spherical bins to cartesian bins
bins_transform = angular_dist_to_euclidean_dist(bins)
bootstraps = []
for i in range(Nbootstraps):
# draw a random sample with N points
ra_R, dec_R = uniform_sphere((min(ra), max(ra)),
(min(dec), max(dec)),
2 * len(ra))
data_R = np.asarray(ra_dec_to_xyz(ra_R, dec_R), order='F').T
if i > 0:
# random sample of the data
ind = np.random.randint(0, data.shape[0], data.shape[0])
data_b = data[ind]
else:
data_b = data
bootstraps.append(two_point(data_b, bins_transform, method=method,
data_R=data_R, random_state=rng))
bootstraps = np.asarray(bootstraps)
corr = np.mean(bootstraps, 0)
corr_err = np.std(bootstraps, 0, ddof=1)
return corr, corr_err, bootstraps
In [ ]:
sklearn_has_two_point
In [ ]:
help(KDTree)
In [ ]:
dataxyz=ra_dec_to_xyz(data['ra'],data['dec'])
In [ ]:
dataxyz=np.asarray(dataxyz)
In [ ]:
dataxyz=dataxyz.transpose()
In [ ]:
dataxyz
In [ ]:
dataxyzR=ra_dec_to_xyz(dataR['ra'],dataR['dec'])
In [ ]:
dataxyzR=np.asarray(dataxyzR)
In [ ]:
dataxyzR=dataxyzR.transpose()
In [ ]:
dataxyzR
In [ ]:
bins=np.arange(0.0,1.05,0.05)
In [ ]:
bins
In [ ]:
In [ ]:
#@pickle_results("tpcf_std.pkl")
tpcf=two_point(dataxyz,bins,method='standard',data_R=dataxyzR, random_state=None)
In [ ]:
tpcf
In [ ]:
plt.plot(bins[1:],tpcf,'ro')
In [ ]:
tpcfam=two_point(dataxyz,bins,method='standard',data_R=None, random_state=None)
In [ ]:
tpcfam
In [ ]:
plt.plot(bins[1:],tpcfam,'bo')
In [ ]:
bins2=np.arange(0.2,0.6,0.02)
In [ ]:
tpcfamb2=two_point(dataxyz,bins2,method='standard',data_R=None, random_state=None)
In [ ]:
plt.plot(bins2[1:],tpcfamb2,'go')
The above doesn't show any BAO feature... It used inbuilt astroML method to generate random catalog... by shuffling the original data's content... That way all of the random points fall in the same survey area and will adhere to all the filtering criteria... the factor or ratio of data pts vs. random pts will be 1... instead of large no. in case if we take existing random catalog or create one
In [ ]:
rng = check_random_state(None)
n_samples, n_features = dataxyz.shape
Nbins = len(bins) - 1
# shuffle all but one axis to get background distribution
data_Rxyz = dataxyz.copy()
print data_Rxyz
for i in range(n_features - 1):
rng.shuffle(data_Rxyz[:, i])
print data_Rxyz
Lets see how it looks with a healpix map
In [ ]:
In [ ]:
NSIDE=512
dr72hpix=hu.HealPix("ring",NSIDE)
In [ ]:
import math as m
def cart2sph(x,y,z):
XsqPlusYsq = x**2 + y**2
r = m.sqrt(XsqPlusYsq + z**2) # r
elev = m.atan2(z,m.sqrt(XsqPlusYsq)) # theta
az = m.atan2(y,x) # phi
return r, elev, az
def cart2sphA(pts):
return np.array([cart2sph(x,y,z) for x,y,z in pts])
def appendSpherical(xyz):
np.hstack((xyz, cart2sphA(xyz)))
In [ ]:
ang=cart2sphA(data_Rxyz)
In [ ]:
ang
In [ ]:
ang.shape
In [ ]:
#ang.resize((105831, 2))
np.squeeze(ang, axis=None)
In [ ]:
help(ang.squeeze)
In [ ]:
ang2=ang[:,1:]
In [ ]:
ang2
In [ ]:
ang2.shape
In [ ]:
ang2[2,0]
In [ ]:
In [ ]:
pixdata = open("./output/pixdatadr72VAGCfullrandam.dat",'w')
pixdata.write("pix \n")
for i in range(0,len(ang2)-1):
#pixdata.write("%f\t" %data['z'][i])
pixdata.write("%d\n" %dr72hpix.ang2pix(ang2[i,0],ang2[i,1]))
pixdata.close()
In [ ]:
pixdata = ascii.read("./output/pixdatadr72VAGCfullrandam.dat")
hpixdata=np.array(np.zeros(hu.nside2npix(NSIDE)))
for j in range(len(pixdata)):
hpixdata[pixdata[j]['pix']]+=1
In [ ]:
hpixdata
In [ ]:
hu.mollview(hpixdata,rot=180)
In [ ]:
hu.orthview(hpixdata)
This method doesnt seem to produce right random catalogs...doing it with ra and dec as follows
In [ ]:
data
In [ ]:
data['z'],data['ra'],data['dec']
In [ ]:
datzradec=np.array([data['z'], data['ra'], data['dec']])
In [ ]:
datzradec
In [ ]:
rng = check_random_state(None)
n_features, n_samples = datzradec.shape
# shuffle all but one axis to get background distribution
data_Rzradec = datzradec.copy()
print data_Rzradec
for i in range(1,n_features):
rng.shuffle(data_Rzradec[:, i])
print data_Rzradec
In [ ]:
min(data_Rzradec[:, 1])
In [ ]:
max(data_Rzradec[:, 1])
In [ ]:
min(data_Rzradec[:, 2])
In [ ]:
max(data_Rzradec[:, 2])
In [ ]:
In [ ]:
min(datzradec[:, 1])
In [ ]:
max(datzradec[:, 1])
In [ ]:
min(datzradec[:, 2])
In [ ]:
max(datzradec[:, 2])
In [ ]:
range(1,3)
In [ ]:
In [ ]:
In [ ]:
help(rng.shuffle)
In [ ]:
n_samples
In [ ]:
n_features
In [ ]:
data_Rzradec
In [ ]:
data_Rzradec[0][2]
In [ ]:
len(data_Rzradec[0][:])
In [ ]:
data_Rzradec[0][:]
In [ ]:
pixdata = open("./output/pixdatadr72VAGCfullrandamrd.dat",'w')
pixdata.write("z\t pix \n")
for i in range(0,len(data_Rzradec[0][:])-1):
pixdata.write("%f\t" %data_Rzradec[0][i])
pixdata.write("%d\n" %dr72hpix.eq2pix(data_Rzradec[1][i],data_Rzradec[2][i]))
pixdata.close()
In [ ]:
pixdata = ascii.read("./output/pixdatadr72VAGCfullrandamrd.dat")
hpixdata=np.array(np.zeros(hu.nside2npix(NSIDE)))
for j in range(len(pixdata)):
hpixdata[pixdata[j]['pix']]+=1
In [ ]:
hpixdata
In [ ]:
hu.mollview(hpixdata,rot=180)
In [ ]:
hu.orthview(hpixdata)
In [ ]:
dataxyz
In [ ]:
dataxyzR1=ra_dec_to_xyz(data_Rzradec[1][:],data_Rzradec[2][:])
In [ ]:
data_Rzradec[1][:]
In [ ]:
In [ ]:
In [ ]:
dataxyzR1
In [ ]:
dataxyzR1=np.asarray(dataxyzR1)
In [ ]:
In [ ]:
dataxyzR1=dataxyzR1.transpose()
In [ ]:
dataxyzR1
In [ ]:
bins=np.arange(0.025,1.025,0.025)
In [ ]:
bins
In [ ]:
In [ ]:
#@pickle_results("tpcf_std.pkl")
tpcf=two_point(dataxyz,bins,method='standard',data_R=dataxyzR1, random_state=None)
In [ ]:
tpcf
In [ ]:
plt.plot(bins[1:],tpcf,'ro')
In [ ]:
bins=np.arange(0.0,1.05,0.05)
In [ ]:
#@pickle_results("tpcf_std.pkl")
tpcf=two_point(dataxyz,bins,method='standard',data_R=dataxyzR1, random_state=None)
In [ ]:
tpcf
In [ ]:
plt.plot(bins[1:],tpcf,'ro')
In [ ]:
btpcf=bootstrap_two_point(dataxyz, bins, Nbootstrap=10,
method='standard', return_bootstraps=False,
random_state=None)
In [ ]:
btpcf
In [ ]:
plt.errorbar(bins[1:],btpcf[0],yerr=btpcf[1],fmt='ro-')
In [ ]:
help(plt.errorbar)
In [ ]:
@pickle_results("tpcf_ls.pkl")
tpcfls=two_point(dataxyz,bins,method='landy-szalay',
data_R=dataxyzR, random_state=None)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
#------------------------------------------------------------
# Set up correlation function computation
# This calculation takes a long time with the bootstrap resampling,
# so we'll save the results.
@pickle_results("correlation_functionsdr72.pkl")
def compute_results(Nbins=16, Nbootstraps=10, method='landy-szalay', rseed=0):
np.random.seed(rseed)
bins = 10 ** np.linspace(np.log10(1. / 60.), np.log10(6), 16)
results = [bins]
for D in [data]:
results += bootstrap_two_point_angular(D['ra'],
D['dec'],
bins=bins,
method=method,
Nbootstraps=Nbootstraps)
return results
(bins, r_corr, r_corr_err, r_bootstraps) = compute_results()
bin_centers = 0.5 * (bins[1:] + bins[:-1])
In [ ]:
bins
In [ ]:
r_corr
In [ ]:
r_corr_err
In [ ]:
r_bootstraps
In [ ]:
#------------------------------------------------------------
# Plot the results
label = '$0.15<z<0.25$\n$N=33813$'
fig = plt.figure(figsize=(6, 6))
plt.xscale('log')
plt.yscale('log')
plt.errorbar(bin_centers, r_corr, r_corr_err,fmt='.k', ecolor='gray', lw=1)
fig.text(0.8, 0.8, label, ha='right', va='top')
plt.xlabel(r'$\theta\ (deg)$')
plt.ylabel(r'$w(\theta)$')
plt.show()
plt.show()
fig.savefig("wth_dr72015025.pdf")
In [ ]:
In [ ]:
data=ascii.read('./input/sdssdr72_sorted_z.dat')
In [ ]:
data
In [ ]:
#m_max = 19
# redshift and magnitude cuts
data = data[data['z'] > 0.05]
data = data[data['z'] < 0.15]
#data = data[data['petroMag_r'] < m_max]
# RA/DEC cuts
RAmin, RAmax = 140, 220
DECmin, DECmax = 5, 45
data = data[data['ra'] < RAmax]
data = data[data['ra'] > RAmin]
data = data[data['dec'] < DECmax]
data = data[data['dec'] > DECmin]
#ur = data['modelMag_u'] - data['modelMag_r']
#flag_red = (ur > 2.22)
#flag_blue = ~flag_red
#datag
print "data size:"
print " total gals: ", len(data)
#print " blue gals:", len(data_blue)
In [ ]:
NSIDE=512
dr72hpix=hu.HealPix("ring",NSIDE)
In [ ]:
pixdata = open("./output/pixdatadr72005015.dat",'w')
pixdata.write("z\t pix \n")
for i in range(0,len(data)-1):
pixdata.write("%f\t" %data['z'][i])
pixdata.write("%d\n" %dr72hpix.eq2pix(data['ra'][i],data['dec'][i]))
pixdata.close()
In [ ]:
pixdata = ascii.read("./output/pixdatadr72005015.dat")
hpixdata=np.array(np.zeros(hu.nside2npix(NSIDE)))
for j in range(len(pixdata)):
hpixdata[pixdata[j]['pix']]+=1
In [ ]:
hpixdata
In [ ]:
hu.mollview(hpixdata,rot=180)
In [ ]:
hu.orthview(hpixdata)
In [ ]:
#------------------------------------------------------------
# Set up correlation function computation
# This calculation takes a long time with the bootstrap resampling,
# so we'll save the results.
@pickle_results("correlation_functionsdr720515.pkl")
def compute_results(Nbins=16, Nbootstraps=10, method='landy-szalay', rseed=0):
np.random.seed(rseed)
bins = 10 ** np.linspace(np.log10(1. / 60.), np.log10(6), 16)
results = [bins]
for D in [data]:
results += bootstrap_two_point_angular(D['ra'],
D['dec'],
bins=bins,
method=method,
Nbootstraps=Nbootstraps)
return results
(bins, r_corr, r_corr_err, r_bootstraps) = compute_results()
bin_centers = 0.5 * (bins[1:] + bins[:-1])
In [ ]:
bins
In [ ]:
r_corr
In [ ]:
r_corr_err
In [ ]:
r_bootstraps
In [ ]:
#------------------------------------------------------------
# Plot the results
label = '$0.05<z<0.15$\n$N=138051$'
fig = plt.figure(figsize=(6, 6))
plt.xscale('log')
plt.yscale('log')
plt.errorbar(bin_centers, r_corr, r_corr_err,fmt='.k', ecolor='gray', lw=1)
fig.text(0.8, 0.8, label, ha='right', va='top')
plt.xlabel(r'$\theta\ (deg)$')
plt.ylabel(r'$w(\theta)$')
plt.show()
fig.savefig("wth_dr72005015.pdf")
In [ ]:
plt.errorbar(bins[0:len(bins)-1],r_corr,r_corr_err)
In [ ]:
data=ascii.read('./input/sdssdr72_sorted_z.dat')
In [ ]:
data
In [ ]:
data['z']
In [ ]:
#m_max = 19
# redshift and magnitude cuts
data = data[data['z'] > 0.05]
data = data[data['z'] <= 0.10]
#data = data[data['petroMag_r'] < m_max]
# RA/DEC cuts
RAmin, RAmax = 140, 220
DECmin, DECmax = 5, 45
data = data[data['ra'] < RAmax]
data = data[data['ra'] > RAmin]
data = data[data['dec'] < DECmax]
data = data[data['dec'] > DECmin]
#ur = data['modelMag_u'] - data['modelMag_r']
#flag_red = (ur > 2.22)
#flag_blue = ~flag_red
#datag
print "data size:"
print " total gals: ", len(data)
#print " blue gals:", len(data_blue)
In [ ]:
NSIDE=512
dr72hpix=hu.HealPix("ring",NSIDE)
In [ ]:
pixdata = open("./output/pixdatadr7200501.dat",'w')
pixdata.write("z\t pix \n")
for i in range(0,len(data)-1):
pixdata.write("%f\t" %data['z'][i])
pixdata.write("%d\n" %dr72hpix.eq2pix(data['ra'][i],data['dec'][i]))
pixdata.close()
In [ ]:
pixdata = ascii.read("./output/pixdatadr7200501.dat")
hpixdata=np.array(np.zeros(hu.nside2npix(NSIDE)))
for j in range(len(pixdata)):
hpixdata[pixdata[j]['pix']]+=1
In [ ]:
hpixdata
In [ ]:
hu.mollview(hpixdata,rot=180)
In [ ]:
hu.orthview(hpixdata)
In [ ]:
#------------------------------------------------------------
# Set up correlation function computation
# This calculation takes a long time with the bootstrap resampling,
# so we'll save the results.
@pickle_results("correlation_functionsdr720501.pkl")
def compute_results(Nbins=16, Nbootstraps=10, method='landy-szalay', rseed=0):
np.random.seed(rseed)
bins = 10 ** np.linspace(np.log10(1. / 60.), np.log10(6), 16)
results = [bins]
for D in [data]:
results += bootstrap_two_point_angular(D['ra'],
D['dec'],
bins=bins,
method=method,
Nbootstraps=Nbootstraps)
return results
(bins, r_corr, r_corr_err, r_bootstraps) = compute_results()
bin_centers = 0.5 * (bins[1:] + bins[:-1])
In [ ]:
bins
In [ ]:
r_corr
In [ ]:
r_corr_err
In [ ]:
r_bootstraps
In [ ]:
#------------------------------------------------------------
# Plot the results
label = '$0.05<z<0.10$\n$N=78939$'
fig = plt.figure(figsize=(6, 6))
plt.xscale('log')
plt.yscale('log')
plt.errorbar(bin_centers, r_corr, r_corr_err,fmt='.k', ecolor='gray', lw=1)
fig.text(0.8, 0.8, label, ha='right', va='top')
plt.xlabel(r'$\theta\ (deg)$')
plt.ylabel(r'$w(\theta)$')
plt.show()
fig.savefig("wth_dr720501.pdf")
In [ ]:
plt.errorbar(bins[0:len(bins)-1],r_corr,r_corr_err)
In [ ]:
In [ ]:
data=ascii.read('./input/sdssdr72_sorted_z.dat')
In [ ]:
data
In [ ]:
data['z']
In [ ]:
#m_max = 19
# redshift and magnitude cuts
data = data[data['z'] > 0.10]
data = data[data['z'] <= 0.15]
#data = data[data['petroMag_r'] < m_max]
# RA/DEC cuts
RAmin, RAmax = 140, 220
DECmin, DECmax = 5, 45
data = data[data['ra'] < RAmax]
data = data[data['ra'] > RAmin]
data = data[data['dec'] < DECmax]
data = data[data['dec'] > DECmin]
#ur = data['modelMag_u'] - data['modelMag_r']
#flag_red = (ur > 2.22)
#flag_blue = ~flag_red
#datag
print "data size:"
print " total gals: ", len(data)
#print " blue gals:", len(data_blue)
In [ ]:
NSIDE=512
dr72hpix=hu.HealPix("ring",NSIDE)
In [ ]:
pixdata = open("./output/pixdatadr72001015.dat",'w')
pixdata.write("z\t pix \n")
for i in range(0,len(data)-1):
pixdata.write("%f\t" %data['z'][i])
pixdata.write("%d\n" %dr72hpix.eq2pix(data['ra'][i],data['dec'][i]))
pixdata.close()
In [ ]:
pixdata = ascii.read("./output/pixdatadr72001015.dat")
hpixdata=np.array(np.zeros(hu.nside2npix(NSIDE)))
for j in range(len(pixdata)):
hpixdata[pixdata[j]['pix']]+=1
In [ ]:
hpixdata
In [ ]:
hu.mollview(hpixdata,rot=180)
In [ ]:
hu.orthview(hpixdata)
In [ ]:
#------------------------------------------------------------
# Set up correlation function computation
# This calculation takes a long time with the bootstrap resampling,
# so we'll save the results.
@pickle_results("correlation_functionsdr72001015.pkl")
def compute_results(Nbins=16, Nbootstraps=10, method='landy-szalay', rseed=0):
np.random.seed(rseed)
bins = 10 ** np.linspace(np.log10(1. / 60.), np.log10(6), 16)
results = [bins]
for D in [data]:
results += bootstrap_two_point_angular(D['ra'],
D['dec'],
bins=bins,
method=method,
Nbootstraps=Nbootstraps)
return results
(bins, r_corr, r_corr_err, r_bootstraps) = compute_results()
bin_centers = 0.5 * (bins[1:] + bins[:-1])
In [ ]:
bins
In [ ]:
r_corr
In [ ]:
r_corr_err
In [ ]:
r_bootstraps
In [ ]:
#------------------------------------------------------------
# Plot the results
label = '$0.10<z<0.15$\n$N=59112$'
fig = plt.figure(figsize=(6, 6))
plt.xscale('log')
plt.yscale('log')
plt.errorbar(bin_centers, r_corr, r_corr_err,fmt='.k', ecolor='gray', lw=1)
fig.text(0.8, 0.8, label, ha='right', va='top')
plt.xlabel(r'$\theta\ (deg)$')
plt.ylabel(r'$w(\theta)$')
plt.show()
fig.savefig("wth_dr7201015.pdf")
In [ ]:
plt.errorbar(bins[0:len(bins)-1],r_corr,r_corr_err)
In [ ]:
hu.mollview(hpixdatab,rot=180)
In [ ]:
hu.orthview(hpixdatab)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
help(hu.mollview)
In [ ]:
from astroML.datasets import fetch_sdss_specgals
from astroML.correlation import bootstrap_two_point_angular
In [ ]:
help(bootstrap_two_point_angular)
In [ ]:
help(astroML.correlation)
In [ ]:
import astroML.correlation
In [ ]:
import sklearn.neighbors
In [ ]:
help(sklearn.neighbors)
Sorted and reduced column set data can now be 'read' to reduce RAM requirements of the table reading.
In [ ]:
sdssdr72=ascii.read('./input/dssdr72_sorted_z.dat')
Create a healpix map with NSIDE=64 (no. of pixels = 49152 as $NPIX=12\times NSIDE^2$) because the no. of galaxies in the survey are less. For higher resolution (later for dr12) we will consider NSIDE=512 or even 1024. For now, we will create a 64 NSIDE map.
In [ ]:
NSIDE=64
dt72hpix=hu.HealPix("ring",NSIDE)
We have data of galaxies with redshifts between 0 and 0.5 ($0<z<0.5$). To look at a time slice/at a certain epoch we need to choose the list of galaxies within a redshift window. As, measurement of redshift has $\pm 0.05$ error. We can bin all the data into redshifts with range limited to 0.05 variation each. So, we have 10 databins with almost identical redshifts. We save each databin in a different file.
In [ ]:
j=0
for i in range(1,17):
pixdata = open("/home/rohin/Desktop/healpix/binned1/pixdata%d_%d.dat"%(NSIDE,i),'w')
pixdata.write("ra\t dec\t z\t pix \n")
#for j in range(len(sdssdr72)):
try:
while sdssdr72[j]['z']<0.03*i:
pixdata.write("%f\t" %sdssdr72[j]['ra'])
pixdata.write("%f\t" %sdssdr72[j]['dec'])
pixdata.write("%f\t" %sdssdr72[j]['z'])
pixdata.write("%d\n" %dt72hpix.eq2pix(sdssdr72[j]['ra'],sdssdr72[j]['dec']))
#print dt72hpix.eq2pix(sdssdr72[j]['ra'],sdssdr72[j]['dec'])
j=j+1
except:
pass
pixdata.close()
In [ ]:
for i in range(1,17):
pixdata = ascii.read("/home/rohin/Desktop/healpix/binned1/pixdata%d_%d.dat"%(NSIDE,i))
mpixdata = open("/home/rohin/Desktop/healpix/binned1/masked/pixdata%d_%d.dat"%(NSIDE,i),'w')
mpixdata.write("ra\t dec\t z\t pix \n")
for j in range((len(pixdata)-1)):
if 100<pixdata[j]['ra']<250:
mpixdata.write("%f\t" %pixdata[j]['ra'])
mpixdata.write("%f\t" %pixdata[j]['dec'])
mpixdata.write("%f\t" %pixdata[j]['z'])
mpixdata.write("%d\n" %pixdata[j]['pix'])
#pixdata.write("/home/rohin/Desktop/healpix/binned1/masked/pixdata_%d.dat"%i,format='ascii')
#print dt72hpix.eq2pix(sdssdr72[j]['ra'],sdssdr72[j]['dec'])
mpixdata.close()
We now, take each databin and assign the total no. of galaxies as the value of each pixel. The following routine will calculate the no. of galaxies by couting the occurence of pixel numbers in the file.
In [ ]:
pixdata = ascii.read("/home/rohin/Desktop/healpix/binned1/masked/pixdata%d_2.dat"%NSIDE)
hpixdata=np.array(np.zeros(hu.nside2npix(NSIDE)))
for j in range(len(pixdata)):
hpixdata[pixdata[j]['pix']]+=1
In [ ]:
hpixdata
In [ ]:
hu.orthview(hpixdata,rot=180)
In [ ]:
pixcl=hu.anafast(hpixdata,lmax=300)
ell = np.arange(len(pixcl))
plt.figure()
plt.plot(ell,np.log(pixcl))
plt.show()
In [ ]:
pixcl=hu.anafast(hpixdata,lmax=300)
ell = np.arange(len(pixcl))
plt.figure()
plt.plot(ell,np.sqrt(ell*(ell+1)*pixcl/(4*math.pi)))
plt.show()
In [ ]:
theta=np.arange(0,np.pi,0.001)
In [ ]:
correldat = np.polynomial.legendre.legval(np.cos(theta),(2*ell+1)*np.absolute(pixcl))/(4*math.pi)
In [ ]:
plt.figure()
plt.plot(theta[0:600]*180/math.pi,correldat[0:600])
plt.show()
In [ ]:
plt.figure()
plt.plot(theta*180/math.pi,correldat)
plt.show()
In [ ]:
randra,randdec=hu.randsphere(2200000)
In [ ]:
randhp=hu.HealPix("RING",NSIDE)
In [ ]:
randhppix=randhp.eq2pix(randra,randdec)
In [ ]:
randpixdat=np.array(np.zeros(hu.nside2npix(NSIDE)))
In [ ]:
for j in range(len(randhppix)):
randpixdat[randhppix[j]]+=1
In [ ]:
randmaphp=hu.mollview(randpixdat)
In [ ]:
randcl=hu.anafast(randpixdat,lmax=300)
ell = np.arange(len(randcl))
plt.figure()
plt.plot(ell,np.sqrt(ell*(ell+1)*randcl/(4*math.pi)))
plt.show()
In [ ]:
correlrand = np.polynomial.legendre.legval(np.cos(theta),(2*ell+1)*np.absolute(randcl))/(4*math.pi)
plt.figure()
plt.plot(theta[0:600]*180/math.pi,correlrand[0:600])
plt.show()
In [ ]:
finalcorrel=correldat-correlrand
plt.figure()
plt.plot(theta[0:600]*180/math.pi,finalcorrel[0:600])
plt.show()
In [ ]:
finalpix=hpixdata-randpixdat
In [ ]:
hu.mollview(finalpix,rot=180)
In [ ]:
cl=hu.anafast(finalpix,lmax=300)
ell = np.arange(len(cl))
plt.figure()
plt.plot(ell,np.sqrt(ell*(ell+1)*cl/(4*math.pi)))
plt.show()
In [ ]:
correlrand = np.polynomial.legendre.legval(np.cos(theta),(2*ell+1)*np.absolute(cl))/(4*math.pi)
plt.figure()
plt.plot(theta[0:600]*180/math.pi,correlrand[0:600])
plt.show()
In [ ]:
finalcl=pixcl-randcl
correlrand = np.polynomial.legendre.legval(np.cos(theta),(2*ell+1)*np.absolute(finalcl))/(4*math.pi)
plt.figure()
plt.plot(theta[0:600]*180/math.pi,correlrand[0:600])
plt.show()
In [ ]:
help(fits)
In [ ]:
data[1].data['z']
In [ ]: