In [1]:
%matplotlib inline
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec
from mpl_toolkits.mplot3d import Axes3D
from astropy.table import Table
from astropy.io import fits as pf
from scipy import interpolate
from scipy.optimize import curve_fit
from scipy import integrate
from scipy import stats
#sys.path.insert(0, '/home/john/densityplot/densityplot')
#from densityplot.hex_scatter import hex_contour as hex_contour
from sklearn.neighbors import KernelDensity as kde
import camb
from camb import model
In [2]:
'''
#Plotting Parameters (Replace with Group code call!)
params = {'legend.fontsize': 16, 'xtick.labelsize': 16, 'ytick.labelsize': 16, 'xtick.major.width':2, 'xtick.minor.width':2, 'ytick.major.width':2, 'ytick.minor.width':2, 'xtick.major.size':8, 'xtick.minor.size':6, 'ytick.major.size':8, 'ytick.minor.size':6}
plt.rcParams.update(params)
plt.rc("axes", linewidth=3.0)
'''
import matplotlib as mpl
def figsze(hscale,
vscale=0.618034,
fig_width_pt = 504.0):
"""Get the fig_width_pt by inserting the textwidth into LaTeX document.
hscale is fraction of text width you want.
vscale is fraction of hscale (defaults to golden ratio)
"""
inches_per_pt = 1.0/72.27 # Convert pt to inch
fig_width = fig_width_pt*inches_per_pt*hscale # width in inches
fig_height = fig_width*vscale # height in inches
fig_size = [fig_width,fig_height]
return fig_size
pgf_with_latex = { # setup matplotlib to use latex for output
"axes.linewidth":1.5, # width of box, 2 is too wide, 1 is too narrow
"pgf.texsystem": "pdflatex", # change this if using xetex or lautex
#"text.usetex": True,
#"text.latex.unicode": True, # use LaTeX to write all text
"font.family": "serif",
"font.serif": [], # blank entries should cause plots to inherit fonts from the document
"font.sans-serif": [],
"font.monospace": [],
"axes.labelsize": 16, # LaTeX default is 10pt font, font size of axis text label
"axes.labelpad" : 6, # Distance between label and axis
"axes.formatter.limits":[-99,99], # use sci notation if log10 of axis range is smaller than first or larger than second.
# GTR: Actually *don't* -- should change the axis label instead. E.g., "Flux Density (10^-17 ergs/s/cm^2)"
# This is a hack b/c there doesn't seem to be an rcParams version of
# axes.ticklabel_format(style='plain')
#"axes.formatter.style":"plain", # Turn off multiplicative offsets (sci notation) to the axes [GTR: Doesn't work]
"axes.formatter.useoffset":False, # Turn off additive offsets to the axes
"font.size": 16,
"legend.fontsize": 12, # Make the legend/label fonts a little smaller
"xtick.labelsize": 16, # Font size of numbers
"ytick.labelsize": 16,
"xtick.direction": "in",
"ytick.direction": "in",
"xtick.minor.visible": True,
"ytick.minor.visible": True,
'xtick.major.width':1,
'xtick.minor.width':1,
'ytick.major.width':1,
'ytick.minor.width':1,
'xtick.major.size':10, # size of tickmarks in points
'xtick.minor.size':5,
'ytick.major.size':10,
'ytick.minor.size':5,
'xtick.major.pad':8, # distance between box and numbers
'ytick.major.pad':8,
"figure.figsize": figsze(1,1), # default fig size of 0.9 textwidth
"pgf.preamble": [
r"\usepackage[utf8x]{inputenc}", # use utf8 fonts because your computer can handle it
r"\usepackage[T1]{fontenc}", # plots will be generated using this preamble
]
}
mpl.rcParams.update(pgf_with_latex)
In [ ]:
In [4]:
#Generate the covariance matrix from the Jackknife file
def make_covariance(Jackfile, warray, RRarray):
JK_xi = []
size = float(len(warray))
jkdat = open(Jackfile,'rw')
#Read the header of the jackknife file and assign the number of jackknife resamples done for the sample
jkhead = jkdat.readline()
print jkhead
jknum= np.float(jkhead.split()[7])
#Pull out the proper info from the file (RR,Xi,etc.) for each jackknife and put into array
RRjk = np.zeros([int(jknum),int(size)])
Xijk = np.zeros([int(jknum),int(size)])
Thjk = np.zeros([int(jknum),int(size)])
num=jkdat.readlines()
row = 0
k=0
for j in range(len(num)):
#At the end of a Jackknife run, I create a new line called 'STEP_DONE' which tells me that my next JK routine is running and separates the two.
if num[j].split()[0] == 'STEP_DONE':
row +=1
k=0
#For each of the JK runs, append the info into arrays
else:
Thjk[row,k] = float(num[j].split()[0])
RRjk[row,k] = float(num[j].split()[3])
Xijk[row,k] = float(num[j].split()[4])
k+=1
C = np.zeros([len(RRjk[0]),len(RRjk[0])])
for m in range(len(RRjk[0])):
for n in range(len(RRjk[0])):
left = np.sqrt(np.asarray(RRjk[:,m])/RRarray[m])
dwa = np.asarray(Xijk[:,m]) - warray[m]
right = np.sqrt(np.asarray(RRjk[:,n])/RRarray[n])
dwb = np.asarray(Xijk[:,n]) - warray[n]
val = left*dwa*right*dwb
C[m,n] = sum(val)
return C
#Generate the covariance matrix from the Jackknife file
def make_JKxi(Jackfile, warray, RRarray):
JK_xi = []
size = float(len(warray))
jkdat = open(Jackfile,'rw')
#Read the header of the jackknife file and assign the number of jackknife resamples done for the sample
jkhead = jkdat.readline()
jknum= np.float(jkhead.split()[7])
#Pull out the proper info from the file (RR,Xi,etc.) for each jackknife and put into array
RRjk = np.zeros([int(jknum),int(size)])
Xijk = np.zeros([int(jknum),int(size)])
Thjk = np.zeros([int(jknum),int(size)])
num=jkdat.readlines()
row = 0
k=0
for j in range(len(num)):
#At the end of a Jackknife run, I create a new line called 'STEP_DONE' which tells me that my next JK routine is running and separates the two.
if num[j].split()[0] == 'STEP_DONE':
row +=1
k=0
#For each of the JK runs, append the info into arrays
else:
Thjk[row,k] = float(num[j].split()[0])
RRjk[row,k] = float(num[j].split()[3])
Xijk[row,k] = float(num[j].split()[4])
k+=1
JK_xi.append(Xijk)
return JK_xi
#Compute chi square
def chisq(dat,mod,cmat):
x = 0
inverse = np.linalg.inv(cmat)
for i in range(len(cmat)):
for j in range(len(cmat)):
df1 = dat[i]-mod[i]
df2 = dat[j]-mod[j]
#icv = cmat[i,j]**-1
icv = inverse[i,j]
x += df1*icv*df2
return x
#Compute the regression matrix
def Regressmat(covar,theta):
R= np.zeros([len(theta),len(theta)])
for i in range(len(theta)):
for j in range(len(theta)):
R[i,j] = covar[i,j]/(np.sqrt(covar[i,i])*np.sqrt(covar[j,j]))
return R
In [5]:
### SPLITTING BY Z Cuting at any imag
#allz= '../Compute_correlation/Extinction_cut_densitycheck15.txt'
allz = '../Compute_correlation2/Final_Clustering_All.txt'
#lowz = '../Compute_correlation/Extinction_cut_densitycheck11.txt'
lowz = '../Compute_correlation2/Final_Clustering_Faint.txt'
#allzJK = '../Compute_correlation/Extinction_cut_densitycheck15_JK.txt'
allzJK = '../Compute_correlation2/Final_Clustering_All_JK.txt'
#lowzJK = '../Compute_correlation/Extinction_cut_densitycheck11_JK.txt'
lowzJK = '../Compute_correlation2/Final_Clustering_Faint_JK.txt'
noextcut = '../Compute_correlation2/Final_Clustering_noextcut.txt'
noextcutJK= '../Compute_correlation2/Final_Clustering_noextcut_JK.txt'
datfiles = [allz,lowz,noextcut]
JKfiles = [allzJK,lowzJK,noextcutJK]
In [6]:
#Loop over the data files list to open and plot the correlation function
separation = []
wnames = []
RRnames = []
DDnames = []
DRnames = []
for infile in datfiles:
fileopen = open(infile,'rw')
header = fileopen.readline()
data = fileopen.readlines()
th = []
w = []
RR = []
DD = []
DR = []
for i in range(len(data)):
t = float(data[i].split()[0]) #Theta
x = float(data[i].split()[4]) #w(theta)
rr = float(data[i].split()[3]) #RR
dd = float(data[i].split()[1]) #DD
dr = float(data[i].split()[2]) #DR
w.append(x)
RR.append(rr)
DD.append(dd)
DR.append(dr)
th.append(t)
#Put the w and RR values into an array to call for Jackknife calculation
wnames.append(w)
RRnames.append(RR)
DDnames.append(DD)
DRnames.append(DR)
separation.append(th)
In [7]:
#Compute Jackknife errors from the counts information
sigma = []
covmat = []
for h in range(len(JKfiles)):
covariance = make_covariance(JKfiles[h], wnames[h], RRnames[h])
sigma.append(np.diag(covariance))
covmat.append(covariance)
#Compute the regression matrix
regmat = []
for h in range(len(wnames)):
regmat.append(Regressmat(covmat[h],th))
print np.shape(JKfiles)
#Get jackknife w's
JK_xi = []
for h in range(len(JKfiles)):
jk = make_JKxi(JKfiles[h], wnames[h], RRnames[h])
JK_xi.append(jk[0])
print np.shape(JK_xi)
print len(separation[0])
In [8]:
#Print values for latex table
pcts = 0
for i in range(len(DDnames[pcts])):
#print separation[0][i], '&', DDnames[0][i], '&', DRnames[0][i], '&', RRnames[0][i], '&', wnames[0][i], '&', sigma[0][i]
print '%.3f & %.1i & %.1i & %.1i & %.3f & %.4f \\\\'%(separation[pcts][i],DDnames[pcts][i],DRnames[pcts][i],RRnames[pcts][i],wnames[pcts][i],sigma[pcts][i]**0.5)
In [9]:
#Following Limber was computed with cosmology (H0,oM,oL)=(70 0.275 0.725)
##This mclimber was done for all objects 2.9<z<5.2 with the final extinction cut
#mclimber_all = [0.0022985455460854311, 0.002190906206060354, 0.0019050009268706861, 0.0014597340543666827, 0.0011930585369182276, 0.0010986709810461117, 0.00094649771865113624, 0.00076558109638027679, 0.00061438375999886519, 0.00044071103670242337, 0.0002922685651020179, 0.0001822212459274252, 9.9271323173412895e-05, 4.6197392904322397e-05, 1.2931319225272703e-05, 4.7453885636605057e-06, -5.7295018738381326e-06, 1.2945112206313357e-06, -2.5558770786635691e-07, 2.0699244496380732e-07]
mclimber_all = [0.0023225429800455637, 0.0021904216380701663, 0.0019131070604180881, 0.0014786904158938835, 0.0011890604062646827, 0.0010940716450842042, 0.00093718770767066959, 0.00076265063034504225, 0.00060422121011192086, 0.00044679016204705126, 0.00030194843029433929, 0.00018166547183924261, 9.4950001759959546e-05, 3.8987292761102685e-05, 1.1656427992710588e-05, 3.8949187589265721e-06, -2.9131135157959959e-06, -1.1865435316490505e-06, -2.3067268951194958e-08, -3.3127561090276543e-07]
#Faint object mclimber with 2.9<z<5.2 with the final extinction cut
#mclimber_low = [0.0023361185250621954, 0.0022064706354088227, 0.0019282439133336759, 0.001482917422530129, 0.0011957113165133413, 0.0011014068340747148, 0.00095946149406393787, 0.0007628377171479422, 0.00061163938639381384, 0.00044805175751709867, 0.00030603189442905508, 0.00018139568437827117, 0.00010018310283104244, 3.6089328262967487e-05, 1.4327521403724255e-05, 2.5924823483705164e-06, -1.7349928197612959e-06, -5.9322324154763784e-07, -1.8647502832907961e-07, -3.2690812395202273e-06]
mclimber_low = [0.0023434681655488164, 0.00220980536106738, 0.0019304367032472322, 0.0014918498380812274, 0.001199407998299017, 0.0011039249516458126, 0.00094488397832085942, 0.00076818403229805582, 0.00060977961027003217, 0.00044978162722766413, 0.00030446958669083705, 0.0001837530609937916, 9.5576443396815135e-05, 3.9697923306540059e-05, 1.1942942460071294e-05, 3.7627680973315105e-06, -2.7558403451147061e-06, -1.3771884816435565e-06, -2.5384149702699815e-07, -1.6631862413132479e-07]
#This mclimber was done for 2.9<z<5.2 with NO extinction cut
mclimber_high = [0.0022431396363423316, 0.0021162407711600136, 0.0018475876015958109, 0.0014284705873674553, 0.0011477142220921229, 0.0010584520197460764, 0.00090569998086366789, 0.00073711419614537932, 0.00058725243468066616, 0.00042998223343848857, 0.00029288227711760491, 0.00017371031323757546, 8.9363275487060279e-05, 3.8516276985481394e-05, 1.1890378726673242e-05, 8.4428691253596875e-07, -4.9585063168790102e-07, -9.6881625471038158e-07, 1.9202865144844596e-07, -1.7224990846038425e-07]
limbera = np.asarray(mclimber_all)*np.pi / 0.7**3
limberl = np.asarray(mclimber_low)*np.pi / 0.7**3
limberh = np.asarray(mclimber_high)*np.pi / 0.7**3
#Interpolate Limber (you get the same function back) to plug in any z in the range (as opposed to set z values)
Limbera = interpolate.interp1d(np.logspace(-1.3,2.5,20),limbera)
Limberl = interpolate.interp1d(np.logspace(-1.3,2.5,20),limberl)
Limberh = interpolate.interp1d(np.logspace(-1.3,2.5,20),limberh)
In [10]:
#Compute different estimators to comapre to LS
EstDP = 103.708998549*np.asarray(DDnames[0])/np.asarray(DRnames[0]) - 1
EstH = np.asarray(DDnames[0]) * np.asarray(RRnames[0])/(np.asarray(DRnames[0]))**2 - 1
print EstDP
print EstH
print wnames[0]
In [11]:
plt.scatter(separation[0],np.asarray(wnames[0])/EstDP,color = 'b', edgecolor = None, s=100, label = r'$\omega_{LS}/\omega_{DP} $')
plt.scatter(separation[0],np.asarray(wnames[0])/EstH,color = 'c', edgecolor = None,s=100, label = r'$\omega_{LS}/\omega_{H} $')
plt.axvline(1,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(30,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.yscale('linear')
plt.ylabel('Estimator Ratio',fontsize = 24)
plt.ylim(-0.5,2)
plt.xscale('log')
plt.xlim(0.1,250)
plt.xlabel(r'$\theta$ (Arcmin)',fontsize = 24)
plt.legend(loc = 3, fontsize = 18, scatterpoints=1)
Out[11]:
In [12]:
#Poisson Errors:
PoiDP = (1.0 + EstDP) / np.sqrt(np.asarray(DDnames[0])/2.0)
PoiDP2 = 103.70899855*np.sqrt(np.asarray(DDnames[0]))/np.asarray(DRnames[0])
PoiH = (1.0 + EstH) / np.sqrt(np.asarray(DDnames[0])/2.0)
PoiLS = (1.0 + np.asarray(wnames[0])) / np.sqrt(np.asarray(DDnames[0])/2.0)
PoiLS1 = (1.0 + np.asarray(wnames[1])) / np.sqrt(np.asarray(DDnames[1])/2.0)
plt.scatter(separation[0],PoiLS/PoiDP,color = 'b', edgecolor = None, s=100, label = r'$\sigma_{LS}/\sigma_{DP} $')
plt.scatter(separation[0],PoiLS/PoiH,color = 'c', edgecolor = None,s=100, label = r'$\sigma_{LS}/\sigma_{H} $')
plt.axvline(1,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(30,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.yscale('linear')
plt.ylabel('Poisson Error Ratio',fontsize = 24)
plt.ylim(-0.5,2)
plt.xscale('log')
plt.xlim(0.1,250)
plt.xlabel(r'$\theta$ (Arcmin)',fontsize = 24)
plt.legend(loc = 3, fontsize = 18, scatterpoints=1)
Out[12]:
In [13]:
#Poisson to JK error
print sigma[0]**0.5
print PoiLS
#plt.scatter(separation[0],sigma[0]**0.5/PoiDP,color = 'b', edgecolor = None, s=100, label = r'$\sigma_{JK}/\sigma_{DP} $')
#plt.scatter(separation[0],sigma[0]**0.5/PoiH,color = 'c', edgecolor = None,s=100, label = r'$\sigma_{JK}/\sigma_{H} $')
plt.scatter(separation[0],sigma[0]**0.5/PoiLS,color = 'g', edgecolor = None,s=100, label = r'$\sigma_{JK}/\sigma_{P} $')
plt.scatter(separation[1],sigma[1]**0.5/PoiLS1,color = 'y', edgecolor = None,s=100, label = r'$\sigma_{JK}/\sigma_{P} $')
plt.axhline(1,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(1,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(30,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.yscale('linear')
plt.ylabel('Error Ratio',fontsize = 24)
plt.ylim(0.5,1.6)
plt.xscale('log')
plt.xlim(0.1,250)
plt.xlabel(r'$\theta$ (Arcmin)',fontsize = 24)
plt.legend(loc = 1, fontsize = 18, scatterpoints=1)
plt.minorticks_on()
#plt.savefig('./Revised_Plots/f15.pdf',bbox_inches='tight')
In [14]:
#New Table 2
pcts = 0
for i in range(len(DDnames[pcts])):
#print separation[0][i], '&', DDnames[0][i], '&', DRnames[0][i], '&', RRnames[0][i], '&', wnames[0][i], '&', sigma[0][i]
print '%.3f & %.1i & %.1i & %.1i & %.3f & %.4f & %.4f \\\\'%(separation[pcts][i],DDnames[pcts][i],DRnames[pcts][i],RRnames[pcts][i],wnames[pcts][i],sigma[pcts][i]**0.5, np.asarray(PoiLS)[i])
In [ ]:
In [15]:
#Replace the Poisson with JK errors
print len(sigma[0])
print len(PoiLS)
maxerr = []
for i in range(len(sigma[0])):
if float(sigma[0][i])**0.5/float(PoiLS[i])>=1.0:
print float(sigma[0][i]), float(PoiLS[i])
maxerr.append(float(sigma[0][i]))
else:
maxerr.append(float(PoiLS[i])**2)
maxerr1 = []
for i in range(len(sigma[1])):
if float(sigma[1][i])**0.5/float(PoiLS1[i])>=1.0:
print float(sigma[1][i]), float(PoiLS1[i])
maxerr1.append(float(sigma[1][i]))
else:
maxerr1.append(float(PoiLS1[i])**2)
sigmatmp = [np.asarray(maxerr),np.asarray(maxerr1),np.asarray(sigma[2])]
In [16]:
#Check to make sure it works!
plt.scatter(separation[0],np.asarray(maxerr)**0.5/PoiLS,color = 'g', edgecolor = None,s=100, label = r'$\sigma_{JK}/\sigma_{H} $')
plt.scatter(separation[1],np.asarray(maxerr1)**0.5/PoiLS1,color = 'y', edgecolor = None,s=100, label = r'$\sigma_{JK}/\sigma_{H} $')
plt.axhline(1,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(1,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(30,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.yscale('linear')
plt.ylabel('Error Ratio',fontsize = 24)
plt.ylim(0.5,1.6)
plt.xscale('log')
plt.xlim(0.1,250)
plt.xlabel(r'$\theta$ (Arcmin)',fontsize = 24)
plt.legend(loc = 1, fontsize = 18, scatterpoints=1)
plt.minorticks_on()
#plt.savefig('./Revised_Plots/f15.pdf',bbox_inches='tight')
In [17]:
#map the changes to the sigma array
sigma = sigmatmp
In [18]:
#Stellar contamination fit using the efficiency of the algorithm (0.86) and the correlation function of stars in the field at large scale (0.09)
def stellar_cont(Limb,b,e):
w_cont = (0.86**2 *b**2 * Limb) + (0.09*(1-0.86)**2) + e
#w_cont = (0.93**2 *b**2 * Limb) + (0.09*(1-0.93)**2) + e
return w_cont
bdx = (np.asarray(separation[0])>1) & (np.asarray(separation[0])<30)& (np.asarray(wnames[0])>0)
#bdx = (np.asarray(separation[0])>1) & (np.asarray(separation[0])<70)& (np.asarray(wnames[0])>-0.1)
allzb, allzcov = curve_fit(stellar_cont, np.asarray(Limbera(separation[0])[bdx]), np.asarray(wnames[0])[bdx], sigma = sigma[0][bdx]**0.5,absolute_sigma=True)#, bounds=([0,0.],[np.inf,np.inf]))
lowzb, lowzcov = curve_fit(stellar_cont, np.asarray(Limberl(separation[1])[bdx]), np.asarray(wnames[1])[bdx], sigma = sigma[1][bdx]**0.5,absolute_sigma=True)#, bounds=([0,0.],[np.inf,np.inf]))
highzb, highzcov = curve_fit(stellar_cont, np.asarray(Limberh(separation[2])[bdx]), np.asarray(wnames[2])[bdx], sigma = sigma[2][bdx]**0.5,absolute_sigma=True)#, bounds=([0,0.],[np.inf,np.inf]))
Las = stellar_cont(Limbera(separation[0]),allzb[0],0)#allzb[1])
Lls = stellar_cont(Limberl(separation[1]),lowzb[0],0)#lowzb[1])
Lhs = stellar_cont(Limberh(separation[2]),highzb[0],0)#highzb[1])
print 'Biases and cross correlations, with contamination:'
print 'allz bias', allzb[0], '+/-', allzcov[0][0]**0.5
print 'allz epsilon',allzb[1], '+/-', allzcov[1][1]**0.5
print ' '
print 'lowz bias',lowzb[0], '+/-', lowzcov[0][0]**0.5
print 'lowz epsilon',lowzb[1], '+/-', lowzcov[1][1]**0.5
stcont_limb = [Las,Lls,Lhs]
In [ ]:
In [ ]:
In [ ]:
In [19]:
#Fit a power law
def powlaw(theta,t0,g):
w_cont = (theta/t0)**(-g)
return w_cont
bdx = (np.asarray(separation[0])>1) & (np.asarray(separation[0])<30)& (np.asarray(wnames[0])>0)
allpow1, apowcov = curve_fit(powlaw, np.asarray(separation[0])[bdx], np.asarray(wnames[0])[bdx], sigma = sigma[0][bdx]**0.5,absolute_sigma=True)
lowpow1, lpowcov = curve_fit(powlaw, np.asarray(separation[1])[bdx], np.asarray(wnames[1])[bdx], sigma = sigma[1][bdx]**0.5,absolute_sigma=True)
highpow1, hpowcov = curve_fit(powlaw, np.asarray(separation[2])[bdx], np.asarray(wnames[2])[bdx], sigma = sigma[2][bdx]**0.5,absolute_sigma=True)
apfit = allpow1
lpfit = lowpow1
hpfit = highpow1
print 'theta0, gamma:'
print ' '
print 'theta =',allpow1[0], '+/-', apowcov[0][0]**0.5
print 'delta =',allpow1[1], '+/-', apowcov[1][1]**0.5
print ' '
print 'theta =',lowpow1[0], '+/-', lpowcov[0][0]**0.5
print 'delta =',lowpow1[1], '+/-', lpowcov[1][1]**0.5
power_idxs = [allpow1,lowpow1,highpow1]
allpow = powlaw(separation[0],allpow1[0],allpow1[1])
lowpow = powlaw(separation[1],lowpow1[0],lowpow1[1])
highpow = powlaw(separation[1],highpow1[0],highpow1[1])
pow_fit = [allpow,lowpow,highpow]
In [20]:
#Chi square of the power law fits
#Need to cast the fitting range to the full matrix
#define a boolean array with the fitting parameters
gdx = (np.asarray(th)>1) & (np.asarray(th)<30)& (np.asarray(wnames[0])>0)
#Change the boolean array to binary; 0=False
int1 = np.asarray(gdx*1)
#Create a binary matrix which reflect where True and False are
gdx2 = int1[:,np.newaxis] * int1
#Convert back to a matrix of booleans
gdx3 = (gdx2 == 1)
chi_pow = []
for h in range(len(wnames)):
chip = chisq(np.asarray(wnames[h])[gdx],pow_fit[h][gdx],covmat[h][gdx3].reshape(len(pow_fit[h][gdx]),len(pow_fit[h][gdx])))
chi_pow.append(chip)
print 'Power law chi^2:', np.asarray(chi_pow)/(len(lowpow[gdx])-2.0), 'for %i DOF'% (len(lowpow[gdx])-2.0)
In [ ]:
In [21]:
#Need to cast the fitting range to the full matrix
#define a boolean array with the fitting parameters
gdx = (np.asarray(th)>1) & (np.asarray(th)<30)& (np.asarray(wnames[0])>0)
#Change the boolean array to binary; 0=False
int1 = np.asarray(gdx*1)
#Create a binary matrix which reflect where True and False are
gdx2 = int1[:,np.newaxis] * int1
#Convert back to a matrix of booleans
gdx3 = (gdx2 == 1)
chi_limb = []
chi_stctm = []
for h in range(len(wnames)):
#chi1 = chisq(np.asarray(wnames[h])[gdx],limb_fit[h][gdx],covmat[h][gdx3].reshape(len(limb_fit[h][gdx]),len(limb_fit[h][gdx])))
chi2 = chisq(np.asarray(wnames[h])[gdx],stcont_limb[h][gdx],covmat[h][gdx3].reshape(len(stcont_limb[h][gdx]),len(stcont_limb[h][gdx])))
#chi_limb.append(chi1)
chi_stctm.append(chi2)
print len(Lls[gdx])
#print 'No stellar contamination chi^2:', np.asarray(chi_limb)/(len(Lls[gdx])-1.0), 'for %i DOF'% (len(Lls[gdx])-1.0)
print 'Stellar contamination chi^2:',np.asarray(chi_stctm)/(len(Lls[gdx])-2.0), 'for %i DOF'% (len(Lls[gdx])-2.0)
In [ ]:
In [22]:
#Integral Constraint
IC = []
for i in range(len(RRnames)):
constraint = np.sum(np.asarray(RRnames[i])*(powlaw(np.asarray(separation[i]),power_idxs[i][0],power_idxs[i][1])))/np.sum(RRnames[i])
IC.append(constraint)
print IC
In [ ]:
In [ ]:
In [41]:
from matplotlib import gridspec
params = {'legend.fontsize': 16, 'xtick.labelsize': 22, 'ytick.labelsize': 22, 'xtick.major.width':2, 'xtick.minor.width':2, 'ytick.major.width':2, 'ytick.minor.width':2, 'xtick.major.size':10, 'xtick.minor.size':6, 'ytick.major.size':10, 'ytick.minor.size':6}
plt.rcParams.update(params)
plt.rc("axes", linewidth=3.0)
ptsz = 150
lwth = 2
#Plot for paper
fig = plt.figure(3,figsize = (10,10))
gs = gridspec.GridSpec(2, 1, height_ratios=[0.75,0.25])
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1],sharex = ax0)
plt.axes(ax0)
plt.scatter(separation[0],wnames[0],s = ptsz,color='#fd8d3c', marker = 'd',edgecolor='None',label=r'QSO Candidates',zorder = 100)
plt.errorbar(separation[0],wnames[0],yerr=sigma[0]**0.5,elinewidth=lwth,fmt=',',color='#fd8d3c')
#Poisson Error
#plt.errorbar(separation[0],wnames[0],yerr=np.sqrt((1.0+np.asarray(wnames[0]))/np.asarray(DDnames[0])),elinewidth=2,fmt=',',color='g')
#No Stellar contamination in fit
#plt.plot(separation[0],La,linestyle = '--', dashes = (8,3,8,3),linewidth = 2,color = '#fd8d3c',label='DM model b=%0.2f'%allzbias[0])
#With Stellar contamination in fit
plt.plot(separation[0],Las,linewidth = lwth,color = '#fd8d3c',label = r'b=%0.2f $\pm$ %0.2f'%(allzb[0],allzcov[0][0]**0.5))
#Power law fit
plt.plot(separation[0],allpow,color='#fd8d3c',linewidth = lwth, linestyle = ':',dashes = (2,2,2,2),label = r'$\theta_0$ =%0.2f ; $\delta$ = %0.2f '%(apfit[0],apfit[1]))
for i in range(len(JK_xi[0])):
if i == len(JK_xi[0])-1:
print 'hi'
plt.plot(separation[0],JK_xi[0][i],color='k',alpha = 0.15, linewidth = 1,label='Jackknife samples')
else:
#plt.scatter(separation[0],JK_xi[0][i],s = 15,color='k',alpha = 0.5, marker = 's',edgecolor='None')
plt.plot(separation[0],JK_xi[0][i],color='k',alpha = 0.15, linewidth = 1)
ax0.set_ylim(5e-4,1.5)
ax0.axvline(1,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
ax0.axvline(30,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
ax0.set_yscale('log')
ax0.set_ylabel(r'$\omega (\theta)$',fontsize = 24)
ax0.legend(fontsize = 18, scatterpoints=1)
plt.axes(ax1)
plt.axhline(1,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(1,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(30,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.scatter(separation[0],np.asarray(wnames[0])/Las,s = ptsz,color='#fd8d3c', marker = 'd',edgecolor='None' )
ax1.yaxis.set_ticks([-2,-1,0,1,2])
ax1.set_ylabel(r'$\omega_m / \omega_{DM}$', fontsize = 20)
plt.subplots_adjust(hspace=.0)
plt.setp(ax0.get_xticklabels(), visible=False)
plt.xlim(0.1,250)
plt.ylim(-2,2)
plt.xscale('log')
plt.xlabel(r'$\theta$ (Arcmin)',fontsize = 24)
#plt.savefig('./Revised_Plots/f9.pdf',bbox_inches='tight')
#plt.savefig('./correlation.png',bbox_inches='tight')
In [24]:
from matplotlib import gridspec
params = {'legend.fontsize': 16, 'xtick.labelsize': 22, 'ytick.labelsize': 22, 'xtick.major.width':2, 'xtick.minor.width':2, 'ytick.major.width':2, 'ytick.minor.width':2, 'xtick.major.size':10, 'xtick.minor.size':6, 'ytick.major.size':10, 'ytick.minor.size':6}
plt.rcParams.update(params)
plt.rc("axes", linewidth=3.0)
ptsz = 200
lwth = 2
#Plot for paper
fig = plt.figure(3,figsize = (10,10))
gs = gridspec.GridSpec(2, 1, height_ratios=[0.75,0.25])
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1],sharex = ax0)
plt.axes(ax0)
plt.scatter(separation[0],wnames[1],s = ptsz,color='#fd8d3c', marker = 'd',edgecolor='None',label=r'Faint QSO Candidates',zorder = 100)
plt.errorbar(separation[0],wnames[1],yerr=sigma[1]**0.5,elinewidth=lwth,fmt=',',color='#fd8d3c')
#Poisson Error
#plt.errorbar(separation[0],wnames[1],yerr=np.sqrt((1.0+np.asarray(wnames[1]))/(np.asarray(DDnames[1])*1.5)),elinewidth=2,fmt=',',color='g')
#No Stellar contamination in fit
#plt.plot(separation[0],La,linestyle = '--', dashes = (8,3,8,3),linewidth = 2,color = '#fd8d3c',label='DM model b=%0.2f'%allzbias[0])
#With Stellar contamination in fit
plt.plot(separation[0],Lls,linewidth = lwth,color = '#fd8d3c',label = r'b=%0.2f $\pm$ %0.2f'%(lowzb[0],lowzcov[0][0]**0.5))
#Power law fit
plt.plot(separation[0],lowpow,color='#fd8d3c',linewidth = lwth, linestyle = ':',dashes = (2,2,2,2),label = r'$\theta_0$ =%0.2f ; $\delta$ = %0.2f '%(lpfit[0],lpfit[1]))
for i in range(len(JK_xi[1])):
if i == len(JK_xi[1])-1:
plt.plot(separation[0],JK_xi[1][i],color='k',alpha = 0.15, linewidth = 1,label='Jackknife samples')
else:
#plt.scatter(separation[0],JK_xi[0][i],s = 15,color='k',alpha = 0.5, marker = 's',edgecolor='None')
plt.plot(separation[0],JK_xi[1][i],color='k',alpha = 0.15, linewidth = 1)
ax0.set_ylim(5e-4,1.5)
ax0.axvline(1,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
ax0.axvline(30,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
ax0.set_yscale('log')
ax0.set_ylabel(r'$\omega (\theta)$',fontsize = 24)
ax0.legend(fontsize = 16, scatterpoints=1)
plt.axes(ax1)
plt.axhline(1,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(1,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(30,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.scatter(separation[1],np.asarray(wnames[1])/Lls,s = ptsz,color='#fd8d3c', marker = 'd',edgecolor='None' )
ax1.yaxis.set_ticks([-2,-1,0,1,2])
ax1.set_ylabel(r'$\omega_m / \omega_{DM}$', fontsize = 20)
plt.subplots_adjust(hspace=.0)
plt.setp(ax0.get_xticklabels(), visible=False)
plt.xlim(0.1,250)
plt.ylim(-2,2)
plt.xscale('log')
plt.xlabel(r'$\theta$ (Arcmin)',fontsize = 24)
#plt.savefig('./Revised_Plots/f10.pdf',bbox_inches='tight')
Out[24]:
In [25]:
from matplotlib import gridspec
params = {'legend.fontsize': 16, 'xtick.labelsize': 22, 'ytick.labelsize': 22, 'xtick.major.width':2, 'xtick.minor.width':2, 'ytick.major.width':2, 'ytick.minor.width':2, 'xtick.major.size':10, 'xtick.minor.size':6, 'ytick.major.size':10, 'ytick.minor.size':6}
plt.rcParams.update(params)
plt.rc("axes", linewidth=3.0)
#Plot for paper
fig = plt.figure(3,figsize = (10,10))
gs = gridspec.GridSpec(2, 1, height_ratios=[0.75,0.25])
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1],sharex = ax0)
plt.axes(ax0)
plt.scatter(separation[0],wnames[0],s = 150,color='#fd8d3c', marker = 'd',edgecolor='None',label=r'Extinction Cut',zorder = 100)
plt.errorbar(separation[0],wnames[0],yerr=sigma[0]**0.5,elinewidth=2,fmt=',',color='#fd8d3c',zorder = 100)
#No Stellar contamination in fit
#plt.plot(separation[0],La,linestyle = '--', dashes = (8,3,8,3),linewidth = 2,color = '#fd8d3c',label='DM model b=%0.2f'%allzbias[0])
#With Stellar contamination in fit
plt.plot(separation[0],Las,linewidth = 2,color = '#fd8d3c',label = r'b=%0.2f $\pm$ %0.2f'%(allzb[0],allzcov[0][0]**0.5))
#Power law fit
plt.plot(separation[0],allpow,color='#fd8d3c',linewidth = 2, linestyle = ':',dashes = (2,2,2,2),label = r'$\theta_0$ =%0.2f ; $\delta$ = %0.2f '%(apfit[0],apfit[1]))
'''
for i in range(len(JK_xi[0])):
if i == len(JK_xi[0])-1:
plt.plot(separation[0],JK_xi[0][i],color='k',alpha = 0.15, linewidth = 1,label='Jackknife samples')
else:
#plt.scatter(separation[0],JK_xi[0][i],s = 15,color='k',alpha = 0.5, marker = 's',edgecolor='None')
plt.plot(separation[0],JK_xi[0][i],color='k',alpha = 0.15, linewidth = 1)
ax0.set_ylim(5e-4,1.5)
'''
plt.scatter(np.asarray(separation[0])+0.1*np.asarray(separation[0]),wnames[2],s = 150,color='g', marker = 'o',edgecolor='None',label=r'No Extinction Cut')
plt.errorbar(np.asarray(separation[0])+0.1*np.asarray(separation[0]),wnames[2],yerr=sigma[2]**0.5,elinewidth=2,fmt=',',color='g')
#No Stellar contamination in fit
#plt.plot(separation[0],La,linestyle = '--', dashes = (8,3,8,3),linewidth = 2,color = '#fd8d3c',label='DM model b=%0.2f'%allzbias[0])
#With Stellar contamination in fit
plt.plot(np.asarray(separation[0])+0.1,Lhs,linewidth = 2,color = 'g',label = r'b=%0.2f $\pm$ %0.2f'%(highzb[0],highzcov[0][0]**0.5))
#Power law fit
plt.plot(np.asarray(separation[0])+0.1,highpow,color='g',linewidth = 2, linestyle = ':',dashes = (2,2,2,2),label = r'$\theta_0$ =%0.2f ; $\delta$ = %0.2f '%(hpfit[0],hpfit[1]))
'''
for i in range(len(JK_xi[0])):
if i == len(JK_xi[0])-1:
plt.plot(separation[0],JK_xi[2][i],color='k',alpha = 0.15, linewidth = 1,label='Jackknife samples')
else:
#plt.scatter(separation[0],JK_xi[0][i],s = 15,color='k',alpha = 0.5, marker = 's',edgecolor='None')
plt.plot(separation[0],JK_xi[2][i],color='k',alpha = 0.15, linewidth = 1)
ax0.set_ylim(5e-4,1.5)
'''
ax0.set_ylim(5e-4,1.5)
ax0.axvline(1,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
ax0.axvline(30,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
ax0.set_yscale('log')
ax0.set_ylabel(r'$\omega (\theta)$',fontsize = 24)
ax0.legend(fontsize = 16, scatterpoints=1)
plt.axes(ax1)
plt.axhline(1,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(1,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(30,linewidth = 2, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.scatter(separation[0],np.asarray(wnames[0])/Las,s = 100,color='#fd8d3c', marker = 'd',edgecolor='None' )
plt.scatter(separation[0]+0.1*np.asarray(separation[0]),np.asarray(wnames[2])/Lhs,s = 100,color='g', marker = 'o',edgecolor='None' )
ax1.yaxis.set_ticks([-2,-1,0,1,2])
ax1.set_ylabel(r'$\omega_m / \omega_{DM}$', fontsize = 20)
plt.subplots_adjust(hspace=.0)
plt.setp(ax0.get_xticklabels(), visible=False)
plt.xlim(0.1,250)
plt.ylim(-2,2)
plt.xscale('log')
plt.xlabel(r'$\theta$ (Arcmin)',fontsize = 24)
#plt.savefig('./Revised_Plots/f13.pdf',bbox_inches='tight')
Out[25]:
In [26]:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure(6,figsize = (10,10))
ax = fig.add_subplot(111, projection='3d')
# Construct arrays for the anchor positions of the 16 bars.
# Note: np.meshgrid gives arrays in (ny, nx) so we use 'F' to flatten xpos,
# ypos in column-major order. For numpy >= 1.7, we could instead call meshgrid
# with indexing='ij'.
X,Y = np.meshgrid(np.log10(th),np.log10(th))
xpos = X.flatten('F')
ypos = Y.flatten('F')
zpos = np.zeros_like(xpos)
# Construct arrays with the dimensions for the 16 bars.
dx = 0.2*np.ones_like(zpos)
dy = dx.copy()
dz = regmat[0].flatten()
ax.bar3d(xpos, ypos, zpos, dx, dy, dz, color='w', zsort='average')
ax.set_zlim(-1,1)
plt.show()
In [27]:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
bdx = (np.asarray(separation[0])>1) & (np.asarray(separation[0])<30)
fig = plt.figure(6,figsize = (10,10))
ax = fig.add_subplot(111, projection='3d')
# Construct arrays for the anchor positions of the 16 bars.
# Note: np.meshgrid gives arrays in (ny, nx) so we use 'F' to flatten xpos,
# ypos in column-major order. For numpy >= 1.7, we could instead call meshgrid
# with indexing='ij'.
X,Y = np.meshgrid(np.log10(np.asarray(separation[0])[bdx]),np.log10(np.asarray(separation[0])[bdx]))
xpos = X.flatten('F')
ypos = Y.flatten('F')
zpos = np.zeros_like(xpos)
# Construct arrays with the dimensions for the 16 bars.
dx = 0.2*np.ones_like(zpos)
dy = dx.copy()
dz = regmat[0][:,bdx][bdx,:].flatten()
ax.bar3d(xpos, ypos, zpos, dx, dy, dz, color='w', zsort='average')
ax.set_zlim(-1,1)
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [28]:
#Biases from different surveys
### He 2017
bhe = 5.93
behe= 1.43
zhe = 3.8
### ROSS 2009
b = [2.06,1.41,1.38,1.45,1.83,2.37,1.92,2.42,2.79,3.62,2.99]
be = [0.03,0.18,0.06,0.38,0.33,0.25,0.5,0.4,0.47,0.49,1.42]
z = [1.27,0.24,0.49,0.80,1.03,1.23,1.41,1.58,1.74,1.92,2.10]
##SHEN 2007 Biases
#sb = [9.8,11.4,13.7]
#sz = [3,3.5,4]
sb = [7.9,14.0]
sbe = [0.8,2.0]
sz = [3.2,4]
### EFTEKHARZADEH
eb = [3.69,3.55,3.57]
ebe = [0.11,0.15,0.09]
ez = [2.297,2.497,2.971]
#### HOPKINS 2007 MODELS
maxdatH07 = open('../Plot_correlation/Hopkins07_clstr_maximal.dat','rw')
defdatH07 = open('../Plot_correlation/Hopkins07_clstr_default.dat','rw')
extdatH07 = open('../Plot_correlation/Hopkins07_clstr_extreme_feedback.dat','rw')
maxH07 = maxdatH07.readlines()
defH07 = defdatH07.readlines()
extH07 = extdatH07.readlines()
zH07 = []
bmH07 = []
bdH07 = []
beH07 = []
b20mH07 = []
b20eH07 = []
b20dH07 = []
for i in range(len(maxH07)):
valm=maxH07[i].split()
vald=defH07[i].split()
vale=extH07[i].split()
zH07.append(np.float(valm[0]))
bmH07.append(np.float(valm[1]))
bdH07.append(np.float(vald[1]))
beH07.append(np.float(vale[1]))
b20mH07.append(np.float(valm[3]))
b20eH07.append(np.float(vale[3]))
b20dH07.append(np.float(vald[3]))
In [42]:
from astropy.cosmology import FlatLambdaCDM
import astropy.units as u
params = {'legend.fontsize': 16, 'xtick.labelsize': 20, 'ytick.labelsize': 20, 'xtick.major.width':3, 'xtick.minor.width':2, 'ytick.major.width':3, 'ytick.minor.width':2, 'xtick.major.size':10, 'xtick.minor.size':5, 'ytick.major.size':10, 'ytick.minor.size':5}
plt.rcParams.update(params)
plt.rc("axes", linewidth=3.0)
#plt.style.use('dark_background')
fig = plt.figure(5,figsize=(10,8))
ax1 = fig.add_subplot(1,1,1)
ax2 = ax1.twiny()
plt.axes(ax1)
plt.scatter(z,b,s=250,c='#0571b0', marker = 'o', edgecolor = 'none', label='Ross 2009 z<2.2')
plt.errorbar(z,b,yerr=be,color='#0571b0', fmt=',',linewidth=3)
plt.scatter(ez,eb,s=250,c='#67a9cf', marker = '^',edgecolor = 'none', label='Eftekharzadeh 2015 2.2<z<2.8')
plt.errorbar(ez,eb,yerr=ebe,color='#67a9cf', fmt=',',linewidth=3)
plt.scatter(sz,sb,s=250,c='#e31a1c', marker = 's',edgecolor = 'none', label='Shen 2007 z>2.9')
plt.errorbar(sz,sb,yerr=sbe,color='#e31a1c', fmt=',',linewidth=3)
plt.scatter(zhe,bhe,s=250,c='m', marker = '^',edgecolor = 'none', label='He 2018 3<z<4')
plt.errorbar(zhe,bhe,yerr=behe,color='m', fmt=',',linewidth=3)
'''
#No stellar contam
plt.scatter(3.479,allzbias, marker = 'd',c='#ca0020', s=100, edgecolor = 'none', label = 'Candidates')
plt.errorbar(3.479,allzbias,yerr=0.22878057**0.5,color='#ca0020', fmt=',',linewidth=3)
plt.scatter(3.143,lowzbias, marker = 'd',c='#ca0020', s=100, edgecolor = 'none')
plt.errorbar(3.143,lowzbias,yerr=0.60656223**0.5,color='#ca0020', fmt=',',linewidth=3)
plt.scatter(3.803,highzbias, marker = 'd',c='#ca0020', s=100, edgecolor = 'none')
plt.errorbar(3.803,highzbias,yerr=0.50286617**0.5,color='#ca0020', fmt=',',linewidth=3)
'''
print allzb[0],allzcov[0][0]
axerror =np.array([[0.34 ,0.24]]).T
ftxerror =np.array([[0.34 ,0.25]]).T
#Stellar contam
plt.scatter(3.38,allzb[0], marker = 'd',c='#fd8d3c', s=350, edgecolor = 'none', label = 'All Candidates')
plt.errorbar(3.38,allzb[0],xerr= axerror, yerr=allzcov[0][0]**0.5,color='#fd8d3c', fmt=',',linewidth=3)
#plt.scatter(3.39,lowzb[0], marker = 'd',c='#a6d96a', s=350, edgecolor = 'none',zorder = 0, label = 'Faint Candidates')
#plt.errorbar(3.39,lowzb[0],xerr=0 ,yerr=lowzcov[0][0]**0.5,color='#a6d96a', fmt=',',linewidth=3,zorder = 0)
#plt.scatter(3.23,highzb[0], marker = 'd',c='#1a9641', s=100, edgecolor = 'none', label = 'Candidates rz match bin')
#plt.errorbar(3.23,highzb[0],xerr=0.0 , yerr=highzcov[0][0]**0.5,color='#1a9641', fmt=',',linewidth=3)
#plt.plot(zH07,b20mH07, linewidth = 2, linestyle='-', color='c',label = "i < 20.2 degeneracy")
plt.plot(zH07,beH07, linewidth = 3, linestyle='-', color='k',label = "Efficient feedback (H07)")
plt.plot(zH07,bmH07, linewidth = 3, linestyle='-.', color='k',dashes = (8,4,2,4,2,4), label = "Maximal growth (H07)")
plt.plot(zH07,bdH07, linewidth = 3, linestyle='--', dashes = (10,5,10,5), color='k',label = "Inefficient feedback (H07)")
#All feedback models
#plt.plot(zH07,b20mH07, linewidth = 2, linestyle=':', c='r',dashes = (3,2,3,2), label = r"All feedback models at $i$<20.2 (H07)")
#plt.plot(zH07,b20eH07, linewidth = 2, linestyle=':', c='r',dashes = (3,2,3,2), label = r"$i$=20.2 (H07)")
#plt.plot(zH07,b20dH07, linewidth = 2, linestyle=':', c='b',dashes = (3,2,3,2), label = r"Feedback models with $i$=20.2 (H07)")
def tick_function(X):
cosmo = FlatLambdaCDM(H0=70 * u.km / u.s / u.Mpc, Tcmb0=2.725 * u.K, Om0=0.274)
ages = cosmo.age(X).value
lab = []
for i in ages:
lab.append(i)
return ["%.2f" % z for z in lab]
ax1.set_xlim(0,4.5)
ax1Ticks = ax1.get_xticks()
ax2Ticks = ax1Ticks
ax2.set_xticks(ax1Ticks)
ax2.set_xbound(ax1.get_xbound())
ax2.set_xticklabels(tick_function(ax2Ticks))
#ax2.set_xticklabels(ax1Ticks)
ax1.set_xlabel('Redshift',fontsize = 24)
ax2.set_xlabel('Age of Universe (Gyr)',fontsize = 24)
plt.ylabel('bias',fontsize = 24)
plt.ylim(0,15)
plt.legend(loc = 2, scatterpoints=1)
plt.minorticks_on()
#plt.savefig('./Revised_Plots/f12.pdf')
#plt.savefig('./Revised_Plots/BZ_all.png')
plt.show()
In [ ]:
In [30]:
#Import SpIES / SHELA data
data = '../Data_Sets/HZLZ_combined_all_hzclassifiers_wphotoz_zspecflg.fits'
obs = pf.open(data)[1].data
Z = obs.zbest
imag = -2.5/np.log(10) * (np.arcsinh((obs.iflux/1e9)/(2*1.8e-10))+np.log(1.8e-10))
gdx = ((Z >= 2.9) & (obs.dec>=-1.2) & (obs.dec<=1.2)& ((obs.ra>=344.4)|(obs.ra<60)))#&(imag >= 20.2))
#gdx = Z>0
#Set up a KDE for dNdz
tmpz = Z[gdx][:, np.newaxis] #change the array from row shape (1) to column shape (1,)
print np.shape(tmpz)
sample_range = np.linspace(min(tmpz[:, 0]), max(tmpz[:, 0]), len(tmpz[:, 0]))[:, np.newaxis]
est = kde(bandwidth=0.1,kernel='epanechnikov') #Set up the Kernel
histkde = est.fit(tmpz).score_samples(sample_range) #fit the kernel to the data and find the density of the grid
#Interpolate (you get the same function back) to plug in any z in the range (as opposed to set z values)
dNdz = interpolate.interp1d(sample_range.flatten(),np.exp(histkde))
print sample_range.flatten()
print 'done'
ZE = np.linspace(min(Z),max(Z),100)
xo=integrate.quad(dNdz,min(sample_range),max(sample_range)) #quad(f(x),xlower,xupper, args)
print xo
print 'median', np.median(Z[gdx])
print 'mean',np.mean(Z[gdx])
print 'max',np.max(Z[gdx])
print 'min',np.min(Z[gdx])
In [31]:
##REDSHIFT SPLITS
In [32]:
#open the angular results from the SHEN 2007 run
shen = './Shen_Angcor.txt'
shenjk = './Shen_Angcor_JK.txt'
datfiles2 = [shen]
JKfiles2 = [shenjk]
#Loop over the data files list to open and plot the correlation function
separation2 = []
wnames2 = []
RRnames2 = []
DDnames2 = []
for infile in datfiles2:
fileopen = open(infile,'rw')
header = fileopen.readline()
data = fileopen.readlines()
th = []
w = []
RR = []
DD = []
for i in range(len(data)):
t = float(data[i].split()[0]) #Theta
x = float(data[i].split()[4]) #w(theta)
rr = float(data[i].split()[3]) #RR
dd = float(data[i].split()[1]) #DD
w.append(x)
RR.append(rr)
DD.append(dd)
th.append(t)
#Put the w and RR values into an array to call for Jackknife calculation
wnames2.append(w)
RRnames2.append(RR)
DDnames2.append(DD)
separation2.append(th)
#Compute Jackknife errors from the counts information
sigma2 = []
covmat2 = []
for h in range(len(JKfiles2)):
covariance = make_covariance(JKfiles[h], wnames[h], RRnames[h])
sigma2.append(np.diag(covariance))
covmat2.append(covariance)
#Compute the regression matrix
regmat2 = []
for h in range(len(wnames2)):
regmat2.append(Regressmat(covmat2[h],th))
print np.shape(regmat2)
#Get jackknife w's
JK_xi2 = []
for h in range(len(JKfiles2)):
jk = make_JKxi(JKfiles2[h], wnames2[h], RRnames2[h])
JK_xi2.append(jk[0])
print np.shape(JK_xi2)
print len(separation2[0])
sth = separation2[0]
sw = wnames2[0]
sjk= sigma2[0]
print len(sth),len(sw),len(sjk)
print DDnames2[0]
In [33]:
#BOSS wtheta which I computed using a subset of Eft 2015
Bth = [0.050323978983, 0.0763011640511, 0.115687744753, 0.175405637024, 0.265949842531, 0.403232871773, 0.611381256445, 0.926975618552, 1.40547945874, 2.13098647839, 3.23099945918, 4.89883798472, 7.42761300449, 11.2617390321, 17.0750368861, 25.8891529834, 39.2531065478, 59.5155188986, 90.2373672122, 136.817801341, 207.443006619]
Bwt = [-0.399949636665, -0.139493890308, -0.446243530879, 0.76311340272, -0.525602284033, -0.573893771779, -0.470477229437, -0.295626322384, 0.207851198524, 0.18424850005, 0.10315515435, 0.0983282856195, 0.0654036946897, 0.0349133960631, 0.0209796601689, 0.0199293887951, 0.014989772678, 0.0125284121898, 0.008932733109, 0.00654631771702, 0.00457786274708]
#The He 2017 points for small scales
heth = [2.05,3.25,5.15,8.15,12.92,20.48,32.46,51.45,81.55,129.24,204.84,324.65,514.53,815.48]
heth = np.asarray(heth)/60.
#Luminous sample
hew = [0,-1,0,-1,-0.33,1.96,-0.85,0.53,0.21,0.15,0.07,-0.07,0.05,-0.005]
hewe= [0,9.76,8.53,1.69,0.86,2.11,0.15,0.32,0.27,0.14,0.08,0.04,0.07,0.03]
#Less Luminous sample
hewll = [0,3.25,2.86,1.58,1.96,-0.11,0.36,-0.002,0.13,0.14,0.06,0.06,0.04,0.02]
hewell = [0,8.74,2.51,1.91,1.79,0.47,0.26,0.18,0.13,0.09,0.04,0.04,0.02,0.02]
In [ ]:
In [34]:
from matplotlib import gridspec
#Plotting Parameters (Replace with Group code call!)
params = {'legend.fontsize': 16, 'xtick.labelsize': 20, 'ytick.labelsize': 20, 'xtick.major.width':3, 'xtick.minor.width':2, 'ytick.major.width':3, 'ytick.minor.width':2, 'xtick.major.size':10, 'xtick.minor.size':5, 'ytick.major.size':10, 'ytick.minor.size':5}
plt.rcParams.update(params)
plt.rc("axes", linewidth=3.0)
#Plot for paper
fig = plt.figure(4,figsize = (8,10))
gs = gridspec.GridSpec(3, 1, height_ratios = [1,1,1], width_ratios=[1])
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1],sharex = ax0)
ax2 = plt.subplot(gs[2],sharex = ax0)
#BOSS
ax0.scatter(Bth,Bwt,s=150,marker = '*', color = 'k',label = r'BOSS $2.2\leq z\leq 2.8$')
#lowz comparison to BOSS
ax0.scatter(separation[0],wnames[0],s = 200,color='#fd8d3c', marker = 'd', edgecolor='None',label='SpIES+SHELA \n'+r'(This work)')
ax0.errorbar(separation[0],wnames[0],yerr=sigma[0]**0.5,elinewidth=1,fmt=',',color='#fd8d3c')
ax0.plot(separation[0],Las,linewidth = 2,color = '#fd8d3c')#,label='Low-z bias=%0.2f'%lowzb[0])
#ax0.plot(separation[1],Ll,linewidth = 2,color = '#a6d96a',label='Low-z bias=%s'%lowzbias)
ax0.plot(separation[0],allpow,color='#fd8d3c',linewidth = 2, linestyle = ':',dashes = (2,2,2,2))#,label = r'$\theta_0$ =%0.2f ; $\delta$ = %0.2f '%(lpfit[0],lpfit[1]))
ax0.set_ylim(5*10**-4,5)
ax0.set_yscale('log')
#ax0.set_ylabel(r'$\omega (\theta)$',fontsize = 18)
ax0.legend(fontsize = 14, scatterpoints=1)
#He less luminous
ax1.scatter(heth,hewll,s = 150,color='c', marker = '^',edgecolor='None',label='HSC $3\leq z\leq 4$'+'\n (CCF)')
ax1.errorbar(heth,hewll,yerr=hewell,elinewidth=1,fmt=',',color='c')
ax1.plot(heth,6.53/60.**0.86*heth**-0.86,color='c',linewidth = 2, linestyle = '--',dashes = (4,2,4,2))
#highz candidates comparison to HE
ax1.scatter(separation[0],wnames[0],s = 200,color='#fd8d3c', marker = 'd', edgecolor='None')
ax1.errorbar(separation[0],wnames[0],yerr=sigma[0]**0.5,elinewidth=1,fmt=',',color='#fd8d3c')
ax1.plot(separation[0],Las,linewidth = 2,color = '#fd8d3c')#,label='Low-z bias=%0.2f'%lowzb[0])
#ax0.plot(separation[1],Ll,linewidth = 2,color = '#a6d96a',label='Low-z bias=%s'%lowzbias)
ax1.plot(separation[0],allpow,color='#fd8d3c',linewidth = 2, linestyle = ':',dashes = (2,2,2,2))#,label = r'$\theta_0$ =%0.2f ; $\delta$ = %0.2f '%(lpfit[0],lpfit[1]))
ax1.set_ylim(5*10**-4,5)
ax1.set_yscale('log')
ax1.set_ylabel(r'$\omega (\theta)$',fontsize = 18)
ax1.legend(fontsize = 14, scatterpoints=1)
#Shen data
ax2.scatter(np.asarray(sth)+0.1*np.asarray(sth),sw,s=150,marker = 'v', color = 'r',label = r'DR5 $2.9\leq z\leq 5.4$ ')
ax2.errorbar(np.asarray(sth)+0.1*np.asarray(sth),sw,yerr=(1.0+np.asarray(sw))/np.sqrt(np.asarray(DDnames2[0])),elinewidth=1,fmt=',',color='r')
#full candidate list
ax2.scatter(separation[0],wnames[0],s = 200,color='#fd8d3c', marker = 'd', edgecolor='None')
ax2.errorbar(separation[0],wnames[0],yerr=sigma[0]**0.5,elinewidth=1,fmt=',',color='#fd8d3c')
ax2.plot(separation[0],Las,linewidth = 2,color = '#fd8d3c')#,label='Low-z bias=%0.2f'%lowzb[0])
#ax0.plot(separation[1],Ll,linewidth = 2,color = '#a6d96a',label='Low-z bias=%s'%lowzbias)
ax2.plot(separation[0],allpow,color='#fd8d3c',linewidth = 2, linestyle = ':',dashes = (2,2,2,2))#,label = r'$\theta_0$ =%0.2f ; $\delta$ = %0.2f '%(lpfit[0],lpfit[1]))
ax2.set_ylim(5*10**-4,5)
ax2.set_yscale('log')
#ax2.set_ylabel(r'$\omega (\theta)$',fontsize = 18)
ax2.legend(fontsize = 14, scatterpoints=1)
plt.subplots_adjust(hspace=.0)
plt.setp(ax0.get_xticklabels(), visible=False)
plt.setp(ax1.get_xticklabels(), visible=False)
plt.xlim(0.1,200)
plt.xscale('log')
plt.xlabel(r'$\theta$ (Arcmin)',fontsize = 18)
#plt.ylabel(r'$\omega (\theta)$',fontsize = 18)
#plt.savefig('./Revised_Plots/f11.pdf',bbox_inches='tight')
#fig.tight_layout()
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [35]:
def stellar_cont(Limb,b,e):
w_cont = (0.86**2 *b**2 * Limb) + (0.09*(1-0.86)**2) + e
#w_cont = (0.93**2 *b**2 * Limb) + (0.09*(1-0.93)**2) + e
return w_cont
def powlaw(theta,t0,g):
w_cont = (theta/t0)**(-g)
return w_cont
In [ ]:
In [ ]:
In [36]:
bdx = (np.asarray(separation[0])>1) & (np.asarray(separation[0])<30)& (np.asarray(wnames[0])>0)
allzb, allzcov = curve_fit(stellar_cont, np.asarray(Limbera(separation[0])[bdx]), np.asarray(wnames[0])[bdx], sigma = sigma[0][bdx]**0.5,absolute_sigma=True)#, bounds=([0,0.],[np.inf,np.inf]))
lowzb, lowzcov = curve_fit(stellar_cont, np.asarray(Limberl(separation[1])[bdx]), np.asarray(wnames[1])[bdx], sigma = sigma[1][bdx]**0.5,absolute_sigma=True)#, bounds=([0,0.],[np.inf,np.inf]))
Las2 = stellar_cont(Limbera(separation[0]),allzb[0],0)#allzb[1])
print allzb[0], '+/-', allzcov[0][0]**0.5
print allzb[1], '+/-', allzcov[1][1]**0.5
print ' '
print lowzb[0], '+/-', lowzcov[0][0]**0.5
print lowzb[1], '+/-', lowzcov[1][1]**0.5
allpow1, apowcov = curve_fit(powlaw, np.asarray(separation[0])[bdx], np.asarray(wnames[0])[bdx], sigma = sigma[0][bdx]**0.5,absolute_sigma=True)
lowpow1, lpowcov = curve_fit(powlaw, np.asarray(separation[1])[bdx], np.asarray(wnames[1])[bdx], sigma = sigma[1][bdx]**0.5,absolute_sigma=True)
allpow = powlaw(separation[0],allpow1[0],allpow1[1])
apfit = allpow1
print ' '
print 'theta =',allpow1[0], '+/-', apowcov[0][0]**0.5
print 'delta =',allpow1[1], '+/-', apowcov[1][1]**0.5
print ' '
print 'theta =',lowpow1[0], '+/-', lpowcov[0][0]**0.5
print 'delta =',lowpow1[1], '+/-', lpowcov[1][1]**0.5
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [37]:
from matplotlib import gridspec
params = {'legend.fontsize': 16, 'xtick.labelsize': 22, 'ytick.labelsize': 22, 'xtick.major.width':2, 'xtick.minor.width':2, 'ytick.major.width':2, 'ytick.minor.width':2, 'xtick.major.size':10, 'xtick.minor.size':6, 'ytick.major.size':10, 'ytick.minor.size':6}
plt.rcParams.update(params)
plt.rc("axes", linewidth=3.0)
ptsz = 300
lwth = 3
#Plot for paper
fig = plt.figure(3,figsize = (10,10))
gs = gridspec.GridSpec(2, 1, height_ratios=[0.75,0.25])
ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1],sharex = ax0)
plt.axes(ax0)
plt.scatter(separation[0],wnames[0],s = ptsz,color='#fd8d3c', marker = 'd',edgecolor='None',label=r'QSO Candidates',zorder = 100)
plt.errorbar(separation[0],wnames[0],yerr=sigma[0]**0.5,elinewidth=lwth,fmt=',',color='#fd8d3c')
#Poisson Error
#plt.errorbar(separation[0],wnames[0],yerr=np.sqrt((1.0+np.asarray(wnames[0]))/np.asarray(DDnames[0])),elinewidth=2,fmt=',',color='g')
#No Stellar contamination in fit
#plt.plot(separation[0],La,linestyle = '--', dashes = (8,3,8,3),linewidth = 2,color = '#fd8d3c',label='DM model b=%0.2f'%allzbias[0])
#With Stellar contamination in fit
plt.plot(separation[0],Las2,linewidth = lwth,color = '#fd8d3c',label = r'b=%0.2f $\pm$ %0.2f'%(allzb[0],allzcov[0][0]**0.5))
#plt.plot(separation[0],stellar_cont(Limbera(separation[0]),5.29,-0.004),linewidth = lwth,color = 'c',label = 'b=5.29')
#Power law fit
plt.plot(separation[0],allpow,color='#fd8d3c',linewidth = lwth, linestyle = ':',dashes = (2,2,2,2),label = r'$\theta_0$ =%0.2f ; $\delta$ = %0.2f '%(apfit[0],apfit[1]))
for i in range(len(JK_xi[0])):
if i == len(JK_xi[0])-1:
plt.plot(separation[0],JK_xi[0][i],color='k',alpha = 0.15, linewidth = 1,label='Jackknife samples')
else:
#plt.scatter(separation[0],JK_xi[0][i],s = 15,color='k',alpha = 0.5, marker = 's',edgecolor='None')
plt.plot(separation[0],JK_xi[0][i],color='k',alpha = 0.15, linewidth = 1)
ax0.set_ylim(5e-4,1.5)
ax0.axvline(1,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
ax0.axvline(30,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
ax0.set_yscale('log')
ax0.set_ylabel(r'$\omega (\theta)$',fontsize = 24)
ax0.legend(fontsize = 18, scatterpoints=1)
plt.axes(ax1)
plt.axhline(1,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(1,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.axvline(30,linewidth = lwth, linestyle = '--', dashes = (8,3,8,3), color = 'k')
plt.scatter(separation[0],np.asarray(wnames[0])/Las2,s = ptsz,color='#fd8d3c', marker = 'd',edgecolor='None' )
ax1.yaxis.set_ticks([-2,-1,0,1,2])
ax1.set_ylabel(r'$\omega_m / \omega_{DM}$', fontsize = 20)
plt.subplots_adjust(hspace=.0)
plt.setp(ax0.get_xticklabels(), visible=False)
plt.xlim(0.1,250)
plt.ylim(-2,2)
plt.xscale('log')
plt.xlabel(r'$\theta$ (Arcmin)',fontsize = 24)
Out[37]:
In [38]:
print allzb[0],allzcov[0][0]
axerror =np.array([[0.34 ,0.24]]).T
ftxerror =np.array([[0.34 ,0.25]]).T
#Stellar contam
plt.scatter(3.38,allzb[0], marker = 'd',c='#fd8d3c', s=200, edgecolor = 'none', label = 'All Candidates')
plt.errorbar(3.38,allzb[0],xerr= axerror, yerr=allzcov[0][0]**0.5,color='#fd8d3c', fmt=',',linewidth=3)
plt.scatter(3.39,lowzb[0], marker = 'd',c='#a6d96a', s=350, edgecolor = 'none',zorder = 0, label = 'Faint Candidates')
plt.errorbar(3.39,lowzb[0],xerr=0 ,yerr=lowzcov[0][0]**0.5,color='#a6d96a', fmt=',',linewidth=3,zorder = 0)
#plt.scatter(3.23,highzb[0], marker = 'd',c='#1a9641', s=100, edgecolor = 'none', label = 'Candidates rz match bin')
#plt.errorbar(3.23,highzb[0],xerr=0.0 , yerr=highzcov[0][0]**0.5,color='#1a9641', fmt=',',linewidth=3)
#plt.plot(zH07,b20mH07, linewidth = 2, linestyle='-', color='c',label = "i < 20.2 degeneracy")
plt.plot(zH07,beH07, linewidth = 3, linestyle='-', color='k',label = "Efficient feedback (H07)")
plt.plot(zH07,bmH07, linewidth = 3, linestyle='-.', color='k',dashes = (8,4,2,4,2,4), label = "Maximal growth (H07)")
plt.plot(zH07,bdH07, linewidth = 3, linestyle='--', dashes = (10,5,10,5), color='k',label = "Inefficient feedback (H07)")
#All feedback models
#plt.plot(zH07,b20mH07, linewidth = 2, linestyle=':', c='r',dashes = (3,2,3,2), label = r"All feedback models at $i$<20.2 (H07)")
#plt.plot(zH07,b20eH07, linewidth = 2, linestyle=':', c='r',dashes = (3,2,3,2), label = r"$i$=20.2 (H07)")
#plt.plot(zH07,b20dH07, linewidth = 2, linestyle=':', c='b',dashes = (3,2,3,2), label = r"Feedback models with $i$=20.2 (H07)")
def tick_function(X):
cosmo = FlatLambdaCDM(H0=70 * u.km / u.s / u.Mpc, Tcmb0=2.725 * u.K, Om0=0.274)
ages = cosmo.age(X).value
lab = []
for i in ages:
lab.append(i)
return ["%.2f" % z for z in lab]
ax1.set_xlim(0,4.5)
ax1Ticks = ax1.get_xticks()
ax2Ticks = ax1Ticks
ax2.set_xticks(ax1Ticks)
ax2.set_xbound(ax1.get_xbound())
ax2.set_xticklabels(tick_function(ax2Ticks))
#ax2.set_xticklabels(ax1Ticks)
ax1.set_xlabel('Redshift',fontsize = 24)
ax2.set_xlabel('Age of Universe (Gyr)',fontsize = 24)
plt.ylabel('bias',fontsize = 24)
plt.ylim(0,15)
plt.xlim(0,4.2)
plt.legend(loc = 2, scatterpoints=1)
plt.minorticks_on()
In [ ]:
In [41]:
#Generating chi square for my different fits
In [39]:
from numpy.linalg import inv
def chi_square1d(data, model, var):
resid = (np.asarray(data) - np.asarray(model))**2
x = sum(resid / np.asarray(var))
return x
def chi_square2d(data, model, covar):
invcov = inv(covar)
resid = data - model
x = sum(sum(np.asarray(resid)[np.newaxis] * invcov * np.asarray(resid)))
return x
def stdev_chisquare(minval, chimat, chival=1):
sdev = np.zeros([np.shape(chimat)[0],np.shape(chimat)[1]])
print np.shape(sdev)
for i in range(len(chimat)):
for j in range(len(chimat[0])):
if chimat[i][j] <= minval + chival:
sdev[i,j] = chimat[i][j]
return sdev
In [40]:
params = {'legend.fontsize': 16, 'xtick.labelsize': 22, 'ytick.labelsize': 22, 'xtick.major.width':2, 'xtick.minor.width':2, 'ytick.major.width':2, 'ytick.minor.width':2, 'xtick.major.size':10, 'xtick.minor.size':6, 'ytick.major.size':10, 'ytick.minor.size':6}
plt.rcParams.update(params)
plt.rc("axes", linewidth=3.0)
In [41]:
#trange = np.arange(0,1.2,0.005)
#grange = np.arange(0.5,2.2,0.005)
trange = np.linspace(0,1.2,400)
grange = np.linspace(0.5,2.2,300)
print len(trange)
chis = np.zeros([len(grange),len(trange)])
bdx = (np.asarray(separation[0])>1) & (np.asarray(separation[0])<30)& (np.asarray(wnames[0])>0)
for i in range(len(trange)):
for j in range(len(grange)):
mod = powlaw(np.asarray(separation[0])[bdx],trange[i],grange[j])
val = chi_square1d(np.asarray(wnames[0])[bdx], mod, np.asarray(sigma[0])[bdx])
chis[j,i] = val
onest = stdev_chisquare(np.amin(chis),chis,1.0)
twost = stdev_chisquare(np.amin(chis),chis,2.71)
X,Y = np.meshgrid(trange,grange)
plt.xlim(np.amin(trange), np.amax(trange))
plt.ylim(np.amin(grange), np.amax(grange))
plt.imshow(chis, extent=(min(trange), max(trange), min(grange), max(grange)),cmap = 'YlOrBr',aspect = 'auto',vmin = 0, vmax = 3,origin = 'lower')#,norm=LogNorm(vmin=0, vmax=30))#, aspect = 'auto')
plt.colorbar(label = r'$\chi^2$')
plt.contour(X[0],Y[:,0],onest,levels=[0],colors = 'k' ,linewidths = 2)
#plt.contour(X[0],Y[:,0],twost,levels=[0],color = 'k', linewidth = 2)
plt.scatter([0.758],[1.447],s=100,color='k',marker = 'x')
plt.xlabel(r'$\theta_0$', fontsize = 18)
plt.ylabel(r'$\delta$', fontsize = 18)
plt.savefig('./Revised_Plots/f16a.pdf',bbox_inches='tight')
plt.show()
In [ ]:
In [42]:
#trange = np.arange(0,1.2,0.005)
#grange = np.arange(0.5,2.2,0.005)
trange = np.linspace(0,1.0,400)
grange = np.linspace(0.4,1.6,300)
print len(trange)
chis = np.zeros([len(grange),len(trange)])
bdx = (np.asarray(separation[1])>1) & (np.asarray(separation[1])<30)& (np.asarray(wnames[1])>0)
for i in range(len(trange)):
for j in range(len(grange)):
mod = powlaw(np.asarray(separation[1])[bdx],trange[i],grange[j])
val = chi_square1d(np.asarray(wnames[1])[bdx], mod, np.asarray(sigma[1])[bdx])
chis[j,i] = val
onest = stdev_chisquare(np.amin(chis),chis,1.0)
X,Y = np.meshgrid(trange,grange)
plt.xlim(np.amin(trange), np.amax(trange))
plt.ylim(np.amin(grange), np.amax(grange))
plt.imshow(chis, extent=(min(trange), max(trange), min(grange), max(grange)),cmap = 'YlOrBr',aspect = 'auto',vmin = 1.5, vmax = 3,origin = 'lower')#,norm=LogNorm(vmin=0, vmax=30))#, aspect = 'auto')
plt.colorbar(label = r'$\chi^2$')
plt.contour(X[0],Y[:,0],onest,levels=[0],colors = 'k' ,linewidths = 2)
plt.scatter([0.373],[0.948],s=100,color='k',marker = 'x')
plt.xlabel(r'$\theta_0$', fontsize = 18)
plt.ylabel(r'$\delta$', fontsize = 18)
plt.savefig('./Revised_Plots/f16b.pdf',bbox_inches='tight')
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [43]:
#Chisquare of the bias function
#brange = np.arange(3,9,0.001)
#erange = np.arange(-0.03,0.03,0.0005)
brange = np.linspace(3,9,300)
erange = np.linspace(-0.03,0.01,600)
print len(erange)
chis = np.zeros([len(erange),len(brange)])
bdx = (np.asarray(separation[0])>1) & (np.asarray(separation[0])<30)& (np.asarray(wnames[0])>0)
for i in range(len(brange)):
for j in range(len(erange)):
mod = stellar_cont(np.asarray(Limbera(separation[0]))[bdx],brange[i],erange[j])
val = chi_square1d(np.asarray(wnames[0])[bdx], mod, np.asarray(sigma[0])[bdx])
chis[j,i] = val
X,Y = np.meshgrid(brange,erange)
onest = stdev_chisquare(np.amin(chis),chis,1.0)
plt.xlim(np.amin(brange), np.amax(brange))
plt.ylim(np.amin(erange), np.amax(erange))
plt.imshow(chis, extent=(min(brange), max(brange), min(erange), max(erange)),cmap = 'YlOrBr',aspect = 'auto',vmin = 1, vmax = 3,origin = 'lower')#,norm=LogNorm(vmin=0, vmax=30))#, aspect = 'auto')
plt.colorbar(label = r'$\chi^2$')
plt.contour(X[0],Y[:,0],onest,levels=[0],colors = 'k' ,linewidths = 2)
plt.scatter([6.488],[-0.0087],s=100,color='k',marker = 'x')
plt.xlabel('Bias', fontsize = 18)
plt.ylabel(r'$\epsilon$', fontsize = 18)
plt.savefig('./Revised_Plots/f17a.pdf',bbox_inches='tight')
plt.show()
In [ ]:
In [44]:
#Chisquare of the bias function for the faint mags
#brange = np.arange(3,9,0.01)
#erange = np.arange(-0.03,0.03,0.0005)
brange = np.linspace(3,9,300)
erange = np.linspace(-0.02,0.03,600)
print len(erange)
chis = np.zeros([len(erange),len(brange)])
bdx = (np.asarray(separation[1])>1) & (np.asarray(separation[1])<30)& (np.asarray(wnames[1])>0)
for i in range(len(brange)):
for j in range(len(erange)):
mod = stellar_cont(np.asarray(Limbera(separation[1]))[bdx],brange[i],erange[j])
val = chi_square1d(np.asarray(wnames[1])[bdx], mod, np.asarray(sigma[1])[bdx])
chis[j,i] = val
onest = stdev_chisquare(np.amin(chis),chis,1.0)
X,Y = np.meshgrid(brange,erange)
plt.xlim(np.amin(brange), np.amax(brange))
plt.ylim(np.amin(erange), np.amax(erange))
plt.imshow(chis, extent=(min(brange), max(brange), min(erange), max(erange)),cmap = 'YlOrBr',aspect = 'auto',vmin = 1.5, vmax = 3,origin = 'lower')#,norm=LogNorm(vmin=0, vmax=30))#, aspect = 'auto')
plt.colorbar(label = r'$\chi^2$')
plt.contour(X[0],Y[:,0],onest,levels=[0],colors = 'k' ,linewidths = 2)
plt.scatter([6.56],[0.0062],s=100,color='k',marker = 'x')
plt.xlabel('Bias', fontsize = 18)
plt.ylabel(r'$\epsilon$', fontsize = 18)
plt.savefig('./Revised_Plots/f17b.pdf',bbox_inches='tight')
plt.show()
In [ ]:
In [ ]: