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 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 [4]:
data=ascii.read("./output/dr12gcmnsrarf.dat")
In [5]:
data
Out[5]:
In [6]:
data.remove_column('z')
data.remove_column('ra')
data.remove_column('dec')
In [7]:
data
Out[7]:
In [8]:
rs=np.array(data['s'])
rrar=np.array(data['rar'])
rdecr=np.array(data['decr'])
In [9]:
dat=np.array([rs,rrar,rdecr])
In [10]:
dat
Out[10]:
In [12]:
dat.reshape(3,len(data))
Out[12]:
In [13]:
dat=dat.transpose()
In [14]:
dat
Out[14]:
In [15]:
# Saving the objects:
with open('dr12gcmnLCDM.pkl', 'w') as f: # Python 3: open(..., 'wb')
pickle.dump(dat, f)
In [3]:
# Getting back the objects:
with open('../pkl/dr12gcmnLCDM.pkl') as f: # Python 3: open(..., 'rb')
dat = pickle.load(f)
dat
Out[3]:
In [4]:
# Getting back the objects:
with open('../pkl/rdr12gcmnLCDM.pkl') as f: # Python 3: open(..., 'rb')
datR = pickle.load(f)
datR
Out[4]:
In [5]:
LCDMmetricsq(dat[0],datR[1])
Out[5]:
In [21]:
%%time
BT_D = BallTree(dat,metric='pyfunc',func=LCDMmetricsq)
with open('BTDdr12gcmnsLCDM.pkl', 'w') as f:
pickle.dump(BT_D,f)
In [8]:
with open('../pkl/BTDdr12gcmnsLCDM.pkl') as f:
BTD = pickle.load(f)
BTD
Out[8]:
In [16]:
bins=np.arange(0.004,0.084,0.004)
bins
Out[16]:
In [17]:
binsq=bins**2
binsq
Out[17]:
In [18]:
%%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/BTDdr12gcmnsDDLCDM.pkl', 'w') as f:
pickle.dump(counts_DD,f)
with open('../pkl/BTDdr12gcmnsDDLCDM.pkl') as f:
counts_DD = pickle.load(f)
counts_DD
In [19]:
DD=np.diff(counts_DD)
DD
Out[19]:
In [ ]:
In [20]:
%%time
start_time=time.time()
counts_DR=BTD.two_point_correlation(datR,binsq)
print counts_DR
end_time=time.time()
tottime=end_time-start_time
print "Total run time:"
print tottime
with open('../pkl/BTDdr12gcmnsDRLCDM.pkl', 'w') as f:
pickle.dump(counts_DR,f)
with open('../pkl/BTDdr12gcmnsDRLCDM.pkl') as f:
counts_DR = pickle.load(f)
counts_DR
In [21]:
counts_DR
Out[21]:
In [22]:
DR=np.diff(counts_DR)
In [23]:
DR
Out[23]:
In [24]:
with open('../pkl/rBTDdr12gcmnsDDLCDM.pkl') as f:
counts_RR = pickle.load(f)
counts_RR
Out[24]:
In [25]:
RR=np.diff(counts_DR)
In [26]:
corrells=(DD+RR-2.0*DR)/RR
In [27]:
plt.plot(bins[1:len(bins)],corrells,'bo-')
Out[27]:
In [29]:
plt.yscale('log')
plt.plot(bins[1:len(bins)],corrells,'bo-')
Out[29]:
In [30]:
plt.yscale('log')
plt.xscale('log')
plt.plot(bins[1:len(bins)],corrells,'bo-')
Out[30]:
In [28]:
corrells
Out[28]:
In [44]:
plt.plot(bins[1:len(bins)],correl,'go-')
Out[44]:
In [45]:
plt.plot(bins[1:len(bins)],bins[1:len(bins)]*bins[1:len(bins)]*correl,'go-')
Out[45]:
In [46]:
plt.plot(bins[2:len(bins)],bins[2:len(bins)]*bins[2:len(bins)]*correl[1:len(bins)],'go-')
Out[46]:
In [47]:
plt.plot(bins[2:len(bins)],correl[1:len(bins)],'go-')
Out[47]:
In [48]:
plt.plot(bins[2:len(bins)],correl[1:len(bins)],'go-')
plt.savefig("correl2x.pdf")
In [53]:
plt.plot(bins[2:len(bins)]*c/1e5,correl[1:len(bins)],'bo-')
plt.savefig("correl2x1.pdf")
In [52]:
plt.yscale('log')
plt.plot(bins[2:len(bins)]*c/1e5,correl[1:len(bins)],'ro-')
Out[52]:
In [54]:
%%time
start_time=time.time()
counts_DR=BTR.two_point_correlation(dat,bins)
print counts_DR
end_time=time.time()
tottime=end_time-start_time
print "Total run time:"
print tottime
with open('BTR200kcDRLCDM.pkl', 'w') as f:
pickle.dump(counts_DR,f)
In [55]:
with open('BTR200kcDRLCDM.pkl') as f:
counts_DR = pickle.load(f)
counts_DR
Out[55]:
In [56]:
DR=np.diff(counts_DR)
In [57]:
corrells=(4.0 * DD - 4.0 * DR + RR) / RR
In [58]:
DR
Out[58]:
In [59]:
corrells
Out[59]:
In [60]:
plt.plot(bins[1:len(bins)],corrells,'go-')
Out[60]:
In [61]:
plt.plot(bins[1:len(bins)],bins[1:len(bins)]*bins[1:len(bins)]*corrells*(c*1e-5)**2,'go-')
Out[61]:
In [72]:
plt.yscale('log')
plt.plot(bins[1:len(bins)]*c/1e5,corrells,'bo-')
plt.savefig("correllsfiglog.pdf")
In [71]:
plt.yscale('log')
plt.plot(bins[2:len(bins)]*c/1e5,corrells[1:len(bins)],'ro-')
plt.savefig("correllslog2x.pdf")
In [73]:
plt.yscale('log')
plt.xscale('log')
plt.plot(bins[1:len(bins)]*c/1e5,corrells,'bo-')
plt.savefig("correllsloglog.pdf")