Skip to content

P
PTAInterpretation
Commit 47c5d4c9 authored by Hannah Middleton's avatar Hannah Middleton
Browse files
Options

some plotting scripts

parent c6a16465
Show whitespace changes
Inline Side-by-side
Showing with 1806 additions and 16 deletions
plotting/GAL_PLOT.py 0 → 100644
View file @ 47c5d4c9
import numpy
import scipy
from numpy.random import choice
import matplotlib.pyplot as plt
import sys
#sys.path.insert(1, '/home/ADF/middlehr/repositories/PTAInference/code/galaxy_model/')
#import mergerrate as mr
sys.path.insert(1, '../models/galaxy_model/')
import mergerratemodel as mr
#import fermif
#import arctanf
green='#009E73' #green
def schechterf(M,M0,phi0,alpha0):
return phi0*numpy.power(M/M0,1.+alpha0)*numpy.exp(-M/M0)
def Mbh(Mstar,alpha,beta):
return alpha+beta*numpy.log10(Mstar/1e11)
def conf_hmz(schains, directory, injection=None, data=None, priorChains=None, modelLoc=None, function=None):
M1 = numpy.linspace(9,12,25)
q = numpy.linspace(0.25,1,10)
z = numpy.linspace(0.01,5.,15)
f = numpy.linspace(-9.,-6.,15)
if data is None:
fbin = 1./(60.*3.154e7) # one over twice tobs in s (half bin width)
else:
fbin = data[0,0]/3.
qdiff = (q[1]-q[0])/2.
zdiff = (z[1]-z[0])/2.
mc = numpy.linspace(5,11,30)
mcdiff = (mc[1]-mc[0])/2.
params = choice(len(schains)-1, 1000, replace=False)
#params = range(len(schains))
h = []
nm = []
nz = []
nsum = []
phig = []
mbh = []
nmerger = []
'''
for k in params:
initpar = dict(Phi0 = schains[k,0], PhiI = schains[k,1], M0 = schains[k,2], alpha = schains[k,3], alphaI = schains[k,4], f0 = schains[k,5], beta = schains[k,7], gamma = schains[k,6], delta = schains[k,8], t0 = schains[k,9], epsilon = schains[k,10], zeta = schains[k,11], eta = schains[k,12], Ms = schains[k,13], theta = schains[k,14], sigma = schains[k,15], e0 = schains[k,16], rho = schains[k,17])
mrate = mr.mergerrate(M1,q,z,f,**initpar)
merger = mrate.hmodelt(fbin=None)
out = numpy.sum(merger[1],axis=1)
#out = numpy.multiply(out,2.*qdiff)
outm = numpy.multiply(out,2.*mrate.zpdiff*mcdiff)
outm = numpy.sum(outm,axis=1)
outz = numpy.sum(out,axis=0)
h.append(merger[0])
nm.append(outm)
nz.append(outz)
nsum.append(numpy.multiply(merger[1],0.55))
#nsum.append(1.1*numpy.sum(outm))
outphi = mrate.output(function='Phi')
phig.append(outphi)
mbh.append(mrate.MBH1)
#nmerger.append(mrate.grid(function=number'))
numpy.savetxt('{0}/h_lines.dat'.format(directory),h,fmt="%f")
numpy.savetxt('{0}/m_lines.dat'.format(directory),nm,fmt="%e")
numpy.savetxt('{0}/z_lines.dat'.format(directory),nz,fmt="%e")
#numpy.savetxt('{0}/nsum.dat'.format(directory),nsum,fmt="%e")
numpy.save('{0}/nsum.npy'.format(directory),nsum)
numpy.save('{0}/phi_g.npy'.format(directory),phig)
numpy.save('{0}/mbh.npy'.format(directory),mbh)
#numpy.save('{0}/number.npy'.format(directory),nmerger)
h = numpy.genfromtxt('{0}/h_lines.dat'.format(directory))
nm = numpy.genfromtxt('{0}/m_lines.dat'.format(directory))
nz = numpy.genfromtxt('{0}/z_lines.dat'.format(directory))
down1,down5,down32,up68,up95,up99,med = numpy.percentile(h,(0.5,5,16,84,95,99.5,50),axis=0)
mdown1,mdown5,mdown32,mup68,mup95,mup99,mmed = numpy.percentile(nm,(0.5,5,16,84,95,99.5,50),axis=0)
zdown1,zdown5,zdown32,zup68,zup95,zup99,zmed = numpy.percentile(nz,(0.5,5,16,84,95,99.5,50),axis=0)
mup99 = numpy.clip(mup99,1e-300,numpy.inf)
mup95 = numpy.clip(mup95,1e-300,numpy.inf)
mup68 = numpy.clip(mup68,1e-300,numpy.inf)
mdown1 = numpy.clip(mdown1,1e-300,numpy.inf)
mdown5 = numpy.clip(mdown5,1e-300,numpy.inf)
mdown32 = numpy.clip(mdown32,1e-300,numpy.inf)
mmed = numpy.clip(mmed,1e-300,numpy.inf)
#sh = numpy.vstack((f,down1,down5,down32,up68,up95,up99,med)).T
#numpy.savetxt("{0}/h.dat".format(directory),sh,fmt='%f')
#sm = numpy.vstack((mc,mdown1,mdown5,mdown32,mup68,mup95,mup99,mmed)).T
#numpy.savetxt("{0}/nm.dat".format(directory),sm,fmt='%e')
#sz = numpy.vstack((z,zdown1,zdown5,zdown32,zup68,zup95,zup99,zmed)).T
#numpy.savetxt("{0}/nz.dat".format(directory),sz,fmt='%f')
'''
# prior plotting
if priorChains is not None:
priorParams = choice(len(priorChains)-1, 1000, replace=False)
#priorParams = range(len(priorChains))
hPrior = []
nmPrior = []
nzPrior = []
#phig_prior = []
#mbh_prior = []
for k in priorParams:
initpar = dict(Phi0 = priorChains[k,0], PhiI = priorChains[k,1], M0 = priorChains[k,2], alpha = priorChains[k,3], alphaI = priorChains[k,4], f0 = priorChains[k,5], beta = priorChains[k,7], gamma = priorChains[k,6], delta = priorChains[k,8], t0 = priorChains[k,9], epsilon = priorChains[k,10], zeta = priorChains[k,11], eta = priorChains[k,12], Ms = priorChains[k,13], theta = priorChains[k,14], sigma = priorChains[k,15], e0 = priorChains[k,16], rho = priorChains[k,17])
mrate = mr.mergerrate(M1,q,z,f,**initpar)
merger = mrate.hmodelt(fbin=None)
out = numpy.sum(merger[1],axis=1)
outm = numpy.multiply(out,mrate.zpdiff/mcdiff)
outm = numpy.sum(outm,axis=1)
outz = numpy.sum(out,axis=0)
hPrior.append(merger[0])
nmPrior.append(outm)
nzPrior.append(outz)
#outphi = mrate.output(function='Phi')
#phig_prior.append(outphi)
#mbh_prior.append(mrate.MBH1)
numpy.savetxt('{0}/h_linesprior.dat'.format(directory),hPrior,fmt="%f")
numpy.savetxt('{0}/m_linesprior.dat'.format(directory),nmPrior,fmt="%e")
numpy.savetxt('{0}/z_linesprior.dat'.format(directory),nzPrior,fmt="%e")
#numpy.save('{0}/phi_g_prior.npy'.format(directory),phig_prior)
#numpy.save('{0}/mbh_prior.npy'.format(directory),mbh_prior)
hPrior = numpy.genfromtxt('{0}/h_linesprior.dat'.format(directory))
nmPrior = numpy.genfromtxt('{0}/m_linesprior.dat'.format(directory))
nzPrior = numpy.genfromtxt('{0}/z_linesprior.dat'.format(directory))
prior1,prior5,prior95,prior99 = numpy.percentile(hPrior,(0.5,5,95,99.5),axis=0)
mprior1,mprior5,mprior95,mprior99 = numpy.percentile(nmPrior,(0.5,5,95,99.5),axis=0)
zprior1,zprior5,zprior95,zprior99 = numpy.percentile(nzPrior,(0.5,5,95,99.5),axis=0)
if injection is not None:
mratei = mr.mergerrate(M1,q,z,f,**injection)
inj0 = mratei.hmodelt(fbin)
inji = mratei.hmodelt()
h0 = inj0[0]
hi = inji[0]
mergeri = numpy.sum(inj0[1],axis=1)
nm0 = numpy.multiply(mergeri,mratei.zpdiff/mcdiff)
nm0 = numpy.sum(nm0,axis=1)
nz0 = numpy.sum(mergeri,axis=0)
f,mc = 10.**f,10.**mc
down1,down5,down32,up68,up95,up99,med = 10.**down1,10.**down5,10.**down32,10.**up68,10.**up95,10.**up99,10.**med
plt.figure(figsize=(8,5))
#plt.locator_params(nbins=4,axis='x')
#plt.locator_params(nbins=6,axis='y')
plt.tick_params(axis='x', labelsize='x-large')
plt.tick_params(axis='y', labelsize='x-large')
if priorChains is not None:
#plt.fill_between(f,10.**prior0,10.**prior100, hatch='.', edgecolor=(0,0,0,1.0),facecolor=(1,1,1,1))
plt.plot(f,10.0**prior1,color=green,linewidth=2,linestyle='dotted')
plt.plot(f,10.0**prior99,color=green,linewidth=2,linestyle='dotted')
plt.plot(f,10.0**prior5,color=green,linewidth=2,linestyle='dashed')
plt.plot(f,10.0**prior95,color=green,linewidth=2,linestyle='dashed')
plt.plot(f,med,color='k',linewidth=2,label='median, 68%, 90%, 95% confidence levels')
#plt.fill_between(f,down1,up99,facecolor='k',alpha=0.25)
plt.fill_between(f,down5,up95,facecolor='k',alpha=0.25)
plt.fill_between(f,down32,up68,facecolor='k',alpha=0.25)
if injection is not None:
plt.plot(f,10.0**h0,color='r',linewidth=2,linestyle='dashdot')
plt.plot(f,10.0**hi,color='r',linewidth=2,linestyle='dotted')
if data is not None:
k = 0
while k < len(data[:,0]):
if data[k,3] < 1.0:
if function == 'fermi':
sigma = fermif.fermis(10.0**data[k,1],0.68)
fermi_up = fermif.fermih(10.0**data[k,1],sigma,0.95)
fermi_down = fermif.fermih(10.0**data[k,1],sigma,0.5)
plt.errorbar(data[k,0],fermi_up,yerr=fermi_up-fermi_down,uplims=True,c='b',linewidth=2)
elif function == 'arctan':
sigma = arctanf.arctans(10.0**data[k,1],0.68)
arctan_up = arctanf.arctanh(10.0**data[k,1],sigma,0.95)
arctan_down = arctanf.arctanh(10.0**data[k,1],sigma,0.5)
plt.errorbar(data[k,0],arctan_up,yerr=arctan_up-arctan_down,uplims=True,c='b',linewidth=2)
elif function == 'ppta':
"""
fraction = 0.8454 #0.8454 is finetuned for specific 4 1.5/T bins ppta datafile
sigma = 10.0**data[k,1]/1.e-15/fraction
fermi_up = 10.0**data[k,1]
fermi_down = sigma*1.e-15
plt.errorbar(data[k,0],fermi_up,yerr=fermi_up-fermi_down,uplims=True,c='b',linewidth=2)
"""
data0 = numpy.genfromtxt('pptaData.dat',skip_footer=1)
plt.plot(data0[:,0],10.0**data0[:,1],c='orangered',linestyle='solid',linewidth=2)
#plt.axvline(6.34e-9,c='orangered',linestyle='dotted',linewidth=2)
plt.plot(6.34e-9,1.e-20*6.34e-9**(-2./3.),c='orangered',marker="*",ms=12)
else: pass
else:
plt.errorbar(data[k,0],10.0**data[k,2],yerr=[[10.0**data[k,2]-10.0**(data[k,2]-data[k,5])],[10.0**(data[k,2]+data[k,5])-10.0**data[k,2]]],c='b',linewidth=2)
k += 1
#plt.plot(data[:,0],10.0**data[:,1],'k',linestyle='dotted',linewidth=2)
plt.errorbar(3.17e-8,1.92e-15,yerr=[[0.55e-15],[0.75e-15]],c='orangered',marker="*",ms=12)
else: pass
#plt.plot(0.,0.,'b',label='simulated data',linewidth=2)
plt.xlabel(r'$ f (\mathrm{Hz})$', fontsize='x-large')
plt.ylabel(r'$ h_c$', fontsize='xx-large')
plt.xlim(1.e-9,1.e-6)
plt.ylim(1.e-18,1.e-13)
plt.xscale('log')
plt.yscale('log')
#plt.legend(loc=3)
plt.savefig("{0}/conf_h.pdf".format(directory),bbox_inches='tight')
plt.savefig("{0}/conf_h.png".format(directory),bbox_inches='tight')
plt.clf()
plt.figure(figsize=(8,5))
#plt.locator_params(nbins=8,axis='x')
#plt.locator_params(nbins=12,axis='y')
plt.tick_params(axis='x', labelsize='x-large')
plt.tick_params(axis='y', labelsize='x-large')
if priorChains is not None:
#plt.fill_between(mc,mprior0,mprior100, hatch='.', edgecolor=(0,0,0,1.0),facecolor=(1,1,1,1))
plt.plot(mc,mprior1,color=green,linewidth=2,linestyle='dotted')
plt.plot(mc,mprior99,color=green,linewidth=2,linestyle='dotted')
plt.plot(mc,mprior5,color=green,linewidth=2,linestyle='dashed')
plt.plot(mc,mprior95,color=green,linewidth=2,linestyle='dashed')
plt.plot(mc,mmed,color='k',linewidth=2,label='median, 68%, 90%, 95% confidence levels')
#plt.fill_between(mc,mdown1,mup99,facecolor='k',alpha=0.25)
plt.fill_between(mc,mdown5,mup95,facecolor='k',alpha=0.25)
plt.fill_between(mc,mdown32,mup68,facecolor='k',alpha=0.25)
if injection is not None:
plt.plot(mc,nm0,color='r',linewidth=2,linestyle='dashdot',label='injected')
plt.xlabel(r'$\mathcal{M} (M_\odot)$', fontsize='x-large')
plt.ylabel(r'$ dN/dVd\log_{10} \mathcal{M} (\mathrm{Mpc}^{-3})$', fontsize='x-large')
plt.xlim(1e6,1e11)
#plt.ylim(1e-9,1e-1)
plt.xscale('log')
plt.yscale('log')
#plt.legend(loc=3)
plt.savefig("{0}/conf_m.pdf".format(directory),bbox_inches='tight')
plt.savefig("{0}/conf_m.png".format(directory),bbox_inches='tight')
if modelLoc is not None:
plotnewmod(modelLoc,'mass')
plt.savefig("{0}/conf_m_model.pdf".format(directory),bbox_inches='tight')
plt.savefig("{0}/conf_m_model.png".format(directory),bbox_inches='tight')
else: pass
plt.clf()
plt.figure(figsize=(8,5))
#plt.locator_params(nbins=8,axis='x')
#plt.locator_params(nbins=12,axis='y')
plt.tick_params(axis='x', labelsize='x-large')
plt.tick_params(axis='y', labelsize='x-large')
if priorChains is not None:
#plt.fill_between(z,zprior0,zprior100, hatch='.', edgecolor=(0,0,0,1.0),facecolor=(1,1,1,1))
plt.plot(z,zprior1,color=green,linewidth=2,linestyle='dotted')
plt.plot(z,zprior99,color=green,linewidth=2,linestyle='dotted')
plt.plot(z,zprior5,color=green,linewidth=2,linestyle='dashed')
plt.plot(z,zprior95,color=green,linewidth=2,linestyle='dashed')
plt.plot(z,zmed,color='k',linewidth=2,label='median, 68%, 90%, 95% confidence levels')
#plt.fill_between(z,zdown1,zup99,facecolor='k',alpha=0.25)
plt.fill_between(z,zdown5,zup95,facecolor='k',alpha=0.25)
plt.fill_between(z,zdown32,zup68,facecolor='k',alpha=0.25)
if injection is not None:
plt.plot(z,nz0,color='r',linewidth=2,linestyle='dashdot',label='injected')
plt.xlabel(r'$z$', fontsize='xx-large')
plt.ylabel(r'$ dN/dVdz (\mathrm{Mpc}^{-3})$', fontsize='x-large')
plt.xlim(0,5)
#plt.ylim(1e-5,1e-1)
plt.yscale('log')
#plt.legend(loc=3)
plt.savefig("{0}/conf_z.pdf".format(directory),bbox_inches='tight')
plt.savefig("{0}/conf_z.png".format(directory),bbox_inches='tight')
if modelLoc is not None:
plotnewmod(modelLoc,'z')
plt.savefig("{0}/conf_z_model.pdf".format(directory),bbox_inches='tight')
plt.savefig("{0}/conf_z_model.png".format(directory),bbox_inches='tight')
else: pass
return
def plotnewmod(dataFile, z_or_m):
newmoddata = numpy.loadtxt('{0}'.format(dataFile,z_or_m))
xvalues=(newmoddata[:,0])
con997_low=(newmoddata[:,1])
con997_high=(newmoddata[:,6])
con2s_low=(newmoddata[:,2])
con2s_high=(newmoddata[:,5])
con1s_low=(newmoddata[:,3])
con1s_high=(newmoddata[:,4])
xnonzero=[]
con997_high_nonzero=[]
con997_low_nonzero=[]
for i in range(len(xvalues)):
if con997_low[i]!=0.0 and con997_high[i]!=0.0:
#if con997_high[i]!=0.0:
xnonzero.append(xvalues[i])
con997_high_nonzero.append(con997_high[i])
con997_low_nonzero.append(con997_low[i])
itop=i
checklow=True
for j in range(len(xvalues)):
if con997_low[j]==0.0 and checklow==True:
con997_low[j]=con997_high[0]
else:
checklow=False
if con997_low[j]==0.0 and checklow==False:
con997_low[j]=con997_high[len(con997_high)-1]
plt.fill_between(xnonzero,con997_low_nonzero,con997_high_nonzero,facecolor='k',edgecolor='k',hatch='XX',alpha=0.5)
return
def conf_g(schains, directory, injection=None, priorChains=None):
M1 = numpy.logspace(9,12,25)
#params = choice(len(schains)-1, 100, replace=False)
params = range(len(schains))
n = []
for k in params:
gsmf = schechterf(M1,10.**schains[k,2],10.**schains[k,0],schains[k,3])
n.append(gsmf)
down1,down5,down32,up68,up95,up99,med = numpy.percentile((n),(0.5,5,16,84,95,99.5,50),axis=0)
up99 = numpy.clip(up99,1e-300,numpy.inf)
up95 = numpy.clip(up95,1e-300,numpy.inf)
up68 = numpy.clip(up68,1e-300,numpy.inf)
down1 = numpy.clip(down1,1e-300,numpy.inf)
down5 = numpy.clip(down5,1e-300,numpy.inf)
down32 = numpy.clip(down32,1e-300,numpy.inf)
med = numpy.clip(med,1e-300,numpy.inf)
s = numpy.vstack((M1,down1,down5,down32,up68,up95,up99,med)).T
numpy.savetxt("{0}/ng.dat".format(directory),s,fmt='%.3e')
# if plotting prior
priorn = []
if priorChains is not None:
priorparams = choice(len(priorChains)-1, 1000, replace=False)
#priorparams = range(len(priorChains))
for k in priorparams:
gsmf = schechterf(M1,10.**priorChains[k,2],10.**priorChains[k,0],priorChains[k,3])
priorn.append(gsmf)
prior0,prior100 = numpy.percentile((priorn),(5.0,95.0),axis=0)
else: pass
if injection is not None:
gsmfi = schechterf(M1,10.**injection['M0'],10.**injection['Phi0'],injection['alpha'])
plt.figure(figsize=(8,5))
#plt.locator_params(nbins=12,axis='x')
#plt.locator_params(nbins=12,axis='y')
plt.tick_params(axis='x', labelsize='x-large')
plt.tick_params(axis='y', labelsize='x-large')
if priorChains is not None:
#plt.fill_between(M1,prior0,prior100, hatch='.', edgecolor=(0,0,0,1.0),facecolor=(1,1,1,1))
plt.plot(M1,prior0,color=green,linewidth=2,linestyle='dashed')
plt.plot(M1,prior100,color=green,linewidth=2,linestyle='dashed')
plt.plot(M1,med,color='k',linewidth=2,label='median, 68%, 90%, 95% confidence levels')
#plt.fill_between(M1,down1,up99,facecolor='k',alpha=0.25)
plt.fill_between(M1,down5,up95,facecolor='k',alpha=0.25)
plt.fill_between(M1,down32,up68,facecolor='k',alpha=0.25)
if injection is not None:
plt.plot(M1,gsmfi,color='r',linewidth=2,linestyle='dashdot',label='injected')
plt.xlabel(r'$M_G (M_\odot)$', fontsize='x-large')
plt.ylabel(r'$\Phi (\mathrm{Mpc}^{-3})$', fontsize='x-large')
plt.xlim(1e9,1e12)
plt.ylim(1e-7,1e-1)
plt.xscale('log')
plt.yscale('log')
#plt.legend(loc=3)
plt.savefig("{0}/conf_g.pdf".format(directory),bbox_inches='tight')
plt.savefig("{0}/conf_g.png".format(directory),bbox_inches='tight')
return
def conf_s(schains, directory, injection=None, priorChains=None):
M1 = numpy.logspace(9,12,25)
#params = choice(len(schains)-1, 100, replace=False)
params = range(len(schains))
n = []
for k in params:
Mb = 10.**Mbh(M1,schains[k,13],schains[k,14])
n.append(Mb)
down1,down5,down32,up68,up95,up99,med = numpy.percentile((n),(0.5,5,16,84,95,99.5,50),axis=0)
s = numpy.vstack((M1,down1,down5,down32,up68,up95,up99,med)).T
numpy.savetxt("{0}/ns.dat".format(directory),s,fmt='%.3e')
# if plotting prior
priorn = []
if priorChains is not None:
priorparams = choice(len(priorChains)-1, 1000, replace=False)
#priorparams = range(len(priorChains))
for k in priorparams:
Mb = 10.**Mbh(M1,priorChains[k,13],priorChains[k,14])
priorn.append(Mb)
prior0,prior100 = numpy.percentile((priorn),(5.0,95.0),axis=0)
else: pass
if injection is not None:
Mbi = 10.**Mbh(M1,injection['Ms'],injection['theta'])
plt.figure(figsize=(8,5))
#plt.locator_params(nbins=12,axis='x')
#plt.locator_params(nbins=12,axis='y')
plt.tick_params(axis='x', labelsize='x-large')
plt.tick_params(axis='y', labelsize='x-large')
if priorChains is not None:
#plt.fill_between(M1,prior0,prior100, hatch='.', edgecolor=(0,0,0,1.0),facecolor=(1,1,1,1))
plt.plot(M1,prior0,color=green,linewidth=2,linestyle='dashed')
plt.plot(M1,prior100,color=green,linewidth=2,linestyle='dashed')
plt.plot(M1,med,color='k',linewidth=2,label='median, 68%, 90%, 95% confidence levels')
#plt.fill_between(M1,down1,up99,facecolor='k',alpha=0.25)
plt.fill_between(M1,down5,up95,facecolor='k',alpha=0.25)
plt.fill_between(M1,down32,up68,facecolor='k',alpha=0.25)
if injection is not None:
plt.plot(M1,Mbi,color='r',linewidth=2,linestyle='dashdot',label='injected')
plt.xlabel(r'$M_{bulge} (M_\odot)$', fontsize='x-large')
plt.ylabel(r'$M_{BH} (M_\odot)$', fontsize='x-large')
plt.xlim(1e9,1e12)
plt.ylim(1e5,1e10)
plt.xscale('log')
plt.yscale('log')
#plt.legend(loc=3)
plt.savefig("{0}/conf_s.pdf".format(directory),bbox_inches='tight')
plt.savefig("{0}/conf_s.png".format(directory),bbox_inches='tight')
return
plotting/agn_corner.py
View file @ 47c5d4c9
......@@ -6,8 +6,57 @@ import matplotlib
matplotlib.rcParams.update({'font.size': 16})
def fiveDCorner(data):
def getAGNPost(nFreqBins):
if nFreqBins==9:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
elif nFreqBins==5:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
elif nFreqBins==3:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
#combine the ns
post0And1 = np.concatenate((post0['logn'], post1['logn']))
logns = np.concatenate((post0And1, post2['logn']))
post0And1 = np.concatenate((post0['alpha'], post1['alpha']))
alphas = np.concatenate((post0And1, post2['alpha']))
post0And1 = np.concatenate((post0['delta'], post1['delta']))
deltas = np.concatenate((post0And1, post2['delta']))
post0And1 = np.concatenate((post0['beta'], post1['beta']))
betas = np.concatenate((post0And1, post2['beta']))
post0And1 = np.concatenate((post0['gamma'], post1['gamma']))
gammas = np.concatenate((post0And1, post2['gamma']))
posteriorCombined = {'logn': logns,
'alpha': alphas,
'delta': deltas,
'beta': betas,
'gamma': gammas}
return posteriorCombined
def fiveDCorner(data,nFBins):
dataMstar = [ np.log10(10.0**(delta+7.0)) for delta in data['delta'] ]
......@@ -23,18 +72,17 @@ def fiveDCorner(data):
figure = corner.corner(d2,labels=parameterLabels,
quantiles=[0.05, 0.5, 0.95])
plt.savefig('corner.png')
plt.show()
plt.savefig('forPaper/agn_fullCorner_fBins{}.pdf'.format(nFBins))
plt.savefig('forPaper/agn_fullCorner_fBins{}.png'.format(nFBins))
return
posterior = np.genfromtxt('posterior.dat', names=True)
fiveDCorner(posterior)
n=9
posterior = getAGNPost(n)#np.genfromtxt('posterior.dat', names=True)
fiveDCorner(posterior,n)
exit()
names = ['logn', 'beta', 'gamma', 'alpha', 'delta']
for i in range(5):
......
plotting/agn_getMlines.py
View file @ 47c5d4c9
......@@ -4,6 +4,55 @@ import astropy.units as u
from astropy.cosmology import Planck18 as cosmo
def getAGNPost(nFreqBins):
if nFreqBins==9:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
elif nFreqBins==5:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
elif nFreqBins==3:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
#combine the ns
post0And1 = np.concatenate((post0['logn'], post1['logn']))
logns = np.concatenate((post0And1, post2['logn']))
post0And1 = np.concatenate((post0['alpha'], post1['alpha']))
alphas = np.concatenate((post0And1, post2['alpha']))
post0And1 = np.concatenate((post0['delta'], post1['delta']))
deltas = np.concatenate((post0And1, post2['delta']))
post0And1 = np.concatenate((post0['beta'], post1['beta']))
betas = np.concatenate((post0And1, post2['beta']))
post0And1 = np.concatenate((post0['gamma'], post1['gamma']))
gammas = np.concatenate((post0And1, post2['gamma']))
posteriorCombined = {'logn': logns,
'alpha': alphas,
'delta': deltas,
'beta': betas,
'gamma': gammas}
saveLinesHere = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{0}/f{0}_plotting_things/'.format(nFreqBins)
return posteriorCombined, saveLinesHere
def zInt(beta,gamma,zlow,zhigh):
......@@ -52,23 +101,25 @@ def dndlog10M(theta,log10Ms):
def main():
pathToRuns = '../../runs/simpleModel/logNormLikeMstar6to10/'
pathToRuns = '../models/agnostic_model/posterior.dat'
nRuns=5
simpleModelData = np.genfromtxt(pathToRuns,names=True)#readPosterior.readPosterior(pathToRuns, nRuns=nRuns)
#pathToRuns = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/nFreq3/posterior.dat'
#pathToRuns = '../models/agnostic_model/posterior.dat'
nFreqBins=9
simpleModelData,resultsDir = getAGNPost(nFreqBins)
# set up M range and write to file
log10MRange = np.linspace(6.0, 11.0, 40)
writeMs = open('tmp2_ms.dat'.format(pathToRuns),'w')
writeMs = open('{}/ms.dat'.format(resultsDir),'w')
for logM in log10MRange:
writeMs.write('{}\n'.format(10.0**logM))
writeMs.close()
# place to write results
writeLines = open('tmp2_mLines.dat'.format(pathToRuns),'w')
writeLines = open('{}/mLines.dat'.format(resultsDir),'w')
# compute lines
for i in range(len(simpleModelData)):
for i in range(len(simpleModelData['logn'])):
thetas = simpleModelData['logn'][i], \
simpleModelData['beta'][i], \
simpleModelData['gamma'][i], \
......
plotting/agn_gethLines.py 0 → 100644
View file @ 47c5d4c9
import numpy as np
from scipy.integrate import quad
import astropy.units as u
from astropy.cosmology import Planck18 as cosmo
import sys
sys.path.insert(1, '../models/agnostic_model/')
import hmodel
def getAGNPost(nFreqBins):
if nFreqBins==9:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
elif nFreqBins==5:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
elif nFreqBins==3:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
#combine the ns
post0And1 = np.concatenate((post0['logn'], post1['logn']))
logns = np.concatenate((post0And1, post2['logn']))
post0And1 = np.concatenate((post0['alpha'], post1['alpha']))
alphas = np.concatenate((post0And1, post2['alpha']))
post0And1 = np.concatenate((post0['delta'], post1['delta']))
deltas = np.concatenate((post0And1, post2['delta']))
post0And1 = np.concatenate((post0['beta'], post1['beta']))
betas = np.concatenate((post0And1, post2['beta']))
post0And1 = np.concatenate((post0['gamma'], post1['gamma']))
gammas = np.concatenate((post0And1, post2['gamma']))
posteriorCombined = {'logn': logns,
'alpha': alphas,
'delta': deltas,
'beta': betas,
'gamma': gammas}
saveLinesHere = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{0}/f{0}_plotting_things/'.format(nFreqBins)
return posteriorCombined, saveLinesHere
def converthc(Ayr,freq):
'''
convert from hc at 1/1yr to whatever frequency we want
'''
# 1/1yr in Hz
freq1yr = 1./(1.*u.year.to(u.s))
# convert
freqRatio = (freq / freq1yr)
hc = Ayr*(freqRatio**(-2./3.))
log10hc = np.log10(float(hc))
return log10hc
posteriors, directory = getAGNPost(9)
nSamples = np.shape(posteriors['delta'])[0]
writeLines = open('./agn_h_nF{}.dat'.format(9),'w')
for i in range(nSamples):
theta = posteriors['logn'][i],\
posteriors['alpha'][i],\
posteriors['delta'][i],\
posteriors['beta'][i],\
posteriors['gamma'][i]
# integration limits
MlowIntLimit = 10.0**6.0
MhighIntLimit = 10.0**11.0
zlowIntLimit = 0.0
zhighIntLimit = 5.0
h1yr = hmodel.hmodel(theta,MlowIntLimit,MhighIntLimit,zlowIntLimit,zhighIntLimit)
logfs = np.linspace(-9.,-6.,15)
hLine = [ converthc(h1yr,10.**logf) for logf in logfs]
for yi in hLine:
writeLines.write('{}\t'.format(yi))
writeLines.write('\n')
plotting/agn_plotChain.py 0 → 100644
View file @ 47c5d4c9
import numpy as np
import matplotlib.pyplot as plt
pathToResults = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/single/'
run='singleChain2'
name = '_s2'
chains = np.genfromtxt(pathToResults+run+'/posterior.dat', names=True)
nparameters = 5
fig, axes = plt.subplots(nparameters, figsize=(10, 10), sharex=True)
labels = (['logn', 'beta', 'gamma', 'alpha', 'delta'])
print(labels)
#labels = ["m", "b", "log(f)"]
for i in range(nparameters):
ax = axes[i]
ax.plot(chains[labels[i]], "k", alpha=0.7)
ax.set_xlim(0, len(chains))
ax.set_ylabel(labels[i],rotation=0)
ax.yaxis.set_label_coords(-0.1, 0.5)
axes[-1].set_xlabel("step number");
plt.savefig('{0}{1}/chains{2}.png'.format(pathToResults,run,name))
plt.show()
plotting/gal_corner.py 0 → 100644
View file @ 47c5d4c9
import corner
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'font.size': 16})
plt.rcParams.update({
"text.usetex": True,
"font.family": "Helvetica"
})
def burninAndThin(samples,burnin,thin):
burnItIn = samples[burnin:]
thinIt = burnItIn[0::thin]
return thinIt
def allCorner(burnin=500,thin=2):
burninP, thinP = 0, 1
nFreqBins = 9
path = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f9/ptmcmc_f{}_1/output/'.format(nFreqBins)
galExt = np.genfromtxt('{}/chain_1.txt'.format(path))
print(galExt)
galParNames = np.loadtxt('/home/ADF/middlehr/repositories/PTAInference/runs/galaxyModel_ext/pars.txt',dtype=str)
prior = np.genfromtxt('/home/ADF/middlehr/repositories/PTAInference/runs/galaxyModel_ext/prior.dat')
burnin = int(0.25*len(galExt))
Phi0Index = list(galParNames).index('Phi0')
Phi0 = burninAndThin(galExt[:,Phi0Index],burnin,thin)
Phi0Prior = burninAndThin(prior[:,Phi0Index],burninP,thinP)
PhiIIndex = list(galParNames).index('PhiI')
PhiI = burninAndThin(galExt[:,PhiIIndex],burnin,thin)
PhiIPrior = burninAndThin(prior[:,PhiIIndex],burninP,thinP)
M0Index = list(galParNames).index('M0')
M0 = burninAndThin(galExt[:,M0Index],burnin,thin)
M0Prior = burninAndThin(prior[:,M0Index],burninP,thinP)
alpha0Index = list(galParNames).index('alpha0')
alpha0 = burninAndThin(galExt[:,alpha0Index],burnin,thin)
alpha0Prior = burninAndThin(prior[:,alpha0Index],burninP,thinP)
alphaIIndex = list(galParNames).index('alphaI')
alphaI = burninAndThin(galExt[:,alphaIIndex],burnin,thin)
alphaIPrior = burninAndThin(prior[:,alphaIIndex],burninP,thinP)
f0Index = list(galParNames).index('f0')
f0 = burninAndThin(galExt[:,f0Index],burnin,thin)
f0Prior = burninAndThin(prior[:,f0Index],burninP,thinP)
betafIndex = list(galParNames).index('betaf')
betaf = burninAndThin(galExt[:,betafIndex],burnin,thin)
betafPrior = burninAndThin(prior[:,betafIndex],burninP,thinP)
alphafIndex = list(galParNames).index('alphaf')
alphaf = burninAndThin(galExt[:,alphafIndex],burnin,thin)
alphafPrior = burninAndThin(prior[:,alphafIndex],burninP,thinP)
gammafIndex = list(galParNames).index('gammaf')
gammaf = burninAndThin(galExt[:,gammafIndex],burnin,thin)
gammafPrior = burninAndThin(prior[:,gammafIndex],burninP,thinP)
tau0Index = list(galParNames).index('tau0')
tau0 = burninAndThin(galExt[:,tau0Index],burnin,thin)
tau0Prior = burninAndThin(prior[:,tau0Index],burninP,thinP)
alphatauIndex = list(galParNames).index('alphatau')
alphatau = burninAndThin(galExt[:,alphatauIndex],burnin,thin)
alphatauPrior = burninAndThin(prior[:,alphatauIndex],burninP,thinP)
betatauIndex = list(galParNames).index('betatau')
betatau = burninAndThin(galExt[:,betatauIndex],burnin,thin)
betatauPrior = burninAndThin(prior[:,betatauIndex],burninP,thinP)
gammatauIndex = list(galParNames).index('gammatau')
gammatau = burninAndThin(galExt[:,gammatauIndex],burnin,thin)
gammatauPrior = burninAndThin(prior[:,gammatauIndex],burninP,thinP)
MstarIndex = list(galParNames).index('Mstar')
Mstar = burninAndThin(galExt[:,MstarIndex],burnin,thin)
MstarPrior = burninAndThin(prior[:,MstarIndex],burninP,thinP)
alphastarIndex = list(galParNames).index('alphastar')
alphastar = burninAndThin(galExt[:,alphastarIndex],burnin,thin)
alphastarPrior = burninAndThin(prior[:,alphastarIndex],burninP,thinP)
epsilonIndex = list(galParNames).index('epsilon')
epsilon = burninAndThin(galExt[:,epsilonIndex],burnin,thin)
epsilonPrior = burninAndThin(prior[:,epsilonIndex],burninP,thinP)
e0Index = list(galParNames).index('e0')
e0 = burninAndThin(galExt[:,e0Index],burnin,thin)
e0Prior = burninAndThin(prior[:,e0Index],burninP,thinP)
zeta0Index = list(galParNames).index('zeta0')
zeta0 = burninAndThin(galExt[:,zeta0Index],burnin,thin)
zeta0Prior = burninAndThin(prior[:,zeta0Index],burninP,thinP)
# do parameter swap
# swap betaf and alphaf
data = np.atleast_1d([Phi0,
PhiI,
M0,
alpha0,
alphaI,
f0,
alphaf,
betaf,
gammaf,
tau0,
alphatau,
betatau,
gammatau,
Mstar,
alphastar,
epsilon,
e0,
zeta0]).T
priorSamples = np.atleast_1d([Phi0Prior,
PhiIPrior,
M0Prior,
alpha0Prior,
alphaIPrior,
f0Prior,
alphafPrior,
betafPrior,
gammafPrior,
tau0Prior,
alphatauPrior,
betatauPrior,
gammatauPrior,
MstarPrior,
alphastarPrior,
epsilonPrior,
e0Prior,
zeta0Prior]).T
parameterLabels = ([r'$\Phi_0$',
r'$\Phi_I$',
r'$\log_{10} M_0$',
r'$\alpha_0$',
r'$\alpha_I$',
r'$f_0$',
r'$\alpha_f$',
r'$\beta_f$',
r'$\gamma_f$',
r'$\tau_0$',
r'$\alpha_{\tau}$',
r'$\beta_{\tau}$',
r'$\gamma_{\tau}$',
r'$\log_{10} M_{*}$',
r'$\alpha_{*}$',
r'$\epsilon$',
r'$e_0$',
r'$\log_{10}\zeta_0$'])
bounds=[(-3.4,-2.4), #Phi0
(-0.6,0.2), #PhiI
(11,11.5), #M0
(-1.5,-1.), # alpha0
(-0.2,0.2), # alphaI
(0.01,0.05), #f0
(-0.5,0.5), #alphaf
(0.,2.), #betaf
(-0.2,0.2), #gammaf
(0.1,5.), #tau0 full range is [0,10]
(-0.5,0.5), #alphatau
(-3.,1.), #betatau
(-0.2,0.2), #gammatau
(7.75,8.75), # M*
(0.9,1.1), #alpha*
(0.2,0.5), #epsilon
(0.01,0.99), #e0
(-2.,2.)] #zeta
figure = corner.corner(data,
levels=(0.5,0.9),
labels=parameterLabels,
quantiles=(0.05,0.50,0.95),
show_titles=True,
title_fmt='.2f',
plot_density=True,
range=bounds,normed=True)
corner.corner(priorSamples,
fig=figure,
no_fill_contours=True,
plot_datapoints=False,
plot_density=False,
levels=([1.]),
color='g',
linewidths=100,
range=bounds,
normed=True)
corner.corner(data,
fig=figure,
levels=(0.5,0.9),
labels=parameterLabels,
quantiles=(0.05,0.50,0.95),
show_titles=True,
title_fmt='.2f',
plot_density=True,
range=bounds,normed=True)
'''
prior = np.genfromtxt('/home/ADF/middlehr/repositories/PTAInference/runs/galaxyModel_ext/prior.dat')
thin=2
mStarBurninPrior = prior[burnin:,mStarIndex]
mStarPrior = mStarBurninPrior[0::thin]
tauBurninPrior = prior[burnin:,tauIndex]
tauPrior = tauBurninPrior[0::thin]
alphaTauBurninPrior = prior[burnin:,alphaTauIndex]
alphaTauPrior = alphaTauBurninPrior[0::thin]
betaTauBurninPrior = prior[burnin:,betaTauIndex]
betaTauPrior = betaTauBurninPrior[0::thin]
priorSamples = np.atleast_1d([tauPrior,alphaTauPrior,betaTauPrior,mStarPrior]).T
print(np.shape(priorSamples))
'''
plt.savefig('forPaper/gal_fullCorner_fBins{}.pdf'.format(nFreqBins))
plt.savefig('forPaper/gal_fullCorner_fBins{}.png'.format(nFreqBins))
plt.show()
return
def selectedCorner(burnin=500,thin=2):
burninP, thinP = 0, 1
nFreqBins = 9
if nFreqBins==9:
path = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f9/ptmcmc_f{}_1/output/'.format(nFreqBins)
elif nFreqBins==5:
path = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f5/ptmcmc_2/output/'.format(nFreqBins)
elif nFreqBins==3:
path = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f3/ptmcmc_f{}_1/output/'.format(nFreqBins)
galExt = np.genfromtxt('{}/chain_1.txt'.format(path))
print(len(galExt))
burnin = int(0.25*len(galExt))
print(burnin)
galParNames = np.loadtxt('/home/ADF/middlehr/repositories/PTAInference/runs/galaxyModel_ext/pars.txt',dtype=str)
prior = np.genfromtxt('/home/ADF/middlehr/repositories/PTAInference/runs/galaxyModel_ext/prior.dat')
print(len(prior))
#tau_0, f_0, M*, e_0, chi, beta_f, alpha_tau, beta_tau
tau0Index = list(galParNames).index('tau0')
tau0 = burninAndThin(galExt[:,tau0Index],burnin,thin)
tau0Prior = burninAndThin(prior[:,tau0Index],burninP,thinP)
f0Index = list(galParNames).index('f0')
f0 = burninAndThin(galExt[:,f0Index],burnin,thin)
f0Prior = burninAndThin(prior[:,f0Index],burninP,thinP)
MstarIndex = list(galParNames).index('Mstar')
Mstar = burninAndThin(galExt[:,MstarIndex],burnin,thin)
MstarPrior = burninAndThin(prior[:,MstarIndex],burninP,thinP)
e0Index = list(galParNames).index('e0')
e0 = burninAndThin(galExt[:,e0Index],burnin,thin)
e0Prior = burninAndThin(prior[:,e0Index],burninP,thinP)
betafIndex = list(galParNames).index('betaf')
betaf = burninAndThin(galExt[:,betafIndex],burnin,thin)
betafPrior = burninAndThin(prior[:,betafIndex],burninP,thinP)
alphatauIndex = list(galParNames).index('alphatau')
alphatau = burninAndThin(galExt[:,alphatauIndex],burnin,thin)
alphatauPrior = burninAndThin(prior[:,alphatauIndex],burninP,thinP)
betatauIndex = list(galParNames).index('betatau')
betatau = burninAndThin(galExt[:,betatauIndex],burnin,thin)
betatauPrior = burninAndThin(prior[:,betatauIndex],burninP,thinP)
M0Index = list(galParNames).index('M0')
M0 = burninAndThin(galExt[:,M0Index],burnin,thin)
M0Prior = burninAndThin(prior[:,M0Index],burninP,thinP)
alpha0Index = list(galParNames).index('alpha0')
alpha0 = burninAndThin(galExt[:,alpha0Index],burnin,thin)
alpha0Prior = burninAndThin(prior[:,alpha0Index],burninP,thinP)
alphaIIndex = list(galParNames).index('alphaI')
alphaI = burninAndThin(galExt[:,alphaIIndex],burnin,thin)
alphaIPrior = burninAndThin(prior[:,alphaIIndex],burninP,thinP)
alphafIndex = list(galParNames).index('alphaf')
alphaf = burninAndThin(galExt[:,alphafIndex],burnin,thin)
alphafPrior = burninAndThin(prior[:,alphafIndex],burninP,thinP)
gammafIndex = list(galParNames).index('gammaf')
gammaf = burninAndThin(galExt[:,gammafIndex],burnin,thin)
gammafPrior = burninAndThin(prior[:,gammafIndex],burninP,thinP)
gammatauIndex = list(galParNames).index('gammatau')
gammatau = burninAndThin(galExt[:,gammatauIndex],burnin,thin)
gammatauPrior = burninAndThin(prior[:,gammatauIndex],burninP,thinP)
alphastarIndex = list(galParNames).index('alphastar')
alphastar = burninAndThin(galExt[:,alphastarIndex],burnin,thin)
alphastarPrior = burninAndThin(prior[:,alphastarIndex],burninP,thinP)
epsilonIndex = list(galParNames).index('epsilon')
epsilon = burninAndThin(galExt[:,epsilonIndex],burnin,thin)
epsilonPrior = burninAndThin(prior[:,epsilonIndex],burninP,thinP)
zeta0Index = list(galParNames).index('zeta0')
zeta0 = burninAndThin(galExt[:,zeta0Index],burnin,thin)
zeta0Prior = burninAndThin(prior[:,zeta0Index],burninP,thinP)
# do parameter swap
# swap betaf and alphaf
data = np.atleast_1d([tau0,
alphatau,
betatau,
f0,
betaf,
Mstar,
e0,
zeta0]).T
'''
M0,
alpha0,
alphaI,
alphaf,
gammaf,
gammatau,
alphastar,
epsilon,
zeta0]).T
'''
priorSamples = np.atleast_1d([tau0Prior,
alphatauPrior,
betatauPrior,
f0Prior,
betafPrior,
MstarPrior,
e0Prior,
zeta0Prior]).T
'''
M0Prior,
alpha0Prior,
alphaIPrior,
alphafPrior,
gammafPrior,
gammatauPrior,
alphastarPrior,
epsilonPrior,
zeta0Prior]).T
'''
parameterLabels = ([r'$\tau_0$',
r'$\alpha_{\tau}$',
r'$\beta_{\tau}$',
r'$f_0$',
r'$\beta_f$',
r'$\log_{10} M_{*}$',
r'$e_0$',
r'$\log_{10}\zeta_0$'])
'''
r'$\log_{10} M_0$',
r'$\alpha_0$',
r'$\alpha_I$',
r'$\alpha_f$',
r'$\gamma_f$',
r'$\gamma_{\tau}$',
r'$\alpha_{*}$',
r'$\epsilon$',
r'$\log_{10}\zeta_0$'])
'''
bounds=[(0.1,5.), #tau0
(-0.5,0.5), #alphatau
(-3.,1.), #betatau
(0.01,0.05), #f0
(0.,2.), #betaf
(7.75,8.75), # M*
(0.01,0.99), #e0
(-2.,2.)] #zeta
figure = corner.corner(data,
levels=(0.5,0.9),
labels=parameterLabels,
quantiles=(0.05,0.50,0.95),
show_titles=True,
title_fmt='.2f',
plot_density=True,
range=bounds)
#smooth1d=True)#,
#normed=True)
corner.corner(priorSamples,
fig=figure,
no_fill_contours=True,
plot_datapoints=False,
plot_density=False,
levels=([1.]),
color='g',
linewidth=100,
range=bounds)
#smooth1d=True)#,
#normed=True)
corner.corner(data,
fig=figure,
levels=(0.5,0.9),
labels=parameterLabels,
quantiles=(0.05,0.50,0.95),
show_titles=True,
title_fmt='.2f',
plot_density=True,
range=bounds)
#smooth1d=True)#,
#normed=True)
'''
prior = np.genfromtxt('/home/ADF/middlehr/repositories/PTAInference/runs/galaxyModel_ext/prior.dat')
thin=2
mStarBurninPrior = prior[burnin:,mStarIndex]
mStarPrior = mStarBurninPrior[0::thin]
tauBurninPrior = prior[burnin:,tauIndex]
tauPrior = tauBurninPrior[0::thin]
alphaTauBurninPrior = prior[burnin:,alphaTauIndex]
alphaTauPrior = alphaTauBurninPrior[0::thin]
betaTauBurninPrior = prior[burnin:,betaTauIndex]
betaTauPrior = betaTauBurninPrior[0::thin]
priorSamples = np.atleast_1d([tauPrior,alphaTauPrior,betaTauPrior,mStarPrior]).T
print(np.shape(priorSamples))
'''
plt.savefig('forPaper/gal_selectedCorner_fBins{}.pdf'.format(nFreqBins))
plt.savefig('forPaper/gal_selectedCorner_fBins{}.png'.format(nFreqBins))
plt.show()
return
allCorner(thin=2)
#selectedCorner(thin=2)
plotting/gal_getMlines.py 0 → 100644
View file @ 47c5d4c9
import numpy as np
import GAL_PLOT
def getchains(nFreqBins):
if nFreqBins==9:
path = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f{0}/ptmcmc_f{0}_1/output/'.format(nFreqBins)
saveResultsHere = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f{0}/f{0}_plotting_things_prior'.format(nFreqBins)
elif nFreqBins==5:
path = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f5/ptmcmc_2/output/'.format(nFreqBins)
elif nFreqBins==3:
path = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f3/ptmcmc_f{}_1/output/'.format(nFreqBins)
readChains = np.genfromtxt('{}/chain_1.txt'.format(path))
burnin = int(0.25*len(readChains))
chains = readChains[burnin:,:]
priorChains = np.genfromtxt('/home/ADF/middlehr/repositories/PTAInference/runs/galaxyModel_ext/prior.dat')
return chains, priorChains, saveResultsHere
chain, prior, directory = getchains(9)
print(chain)
print(prior)
GAL_PLOT.conf_hmz(chain,directory,priorChains=prior)
plotting/gal_hlines.py 0 → 100644
View file @ 47c5d4c9
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'font.size': 16})
plt.rcParams.update({
"text.usetex": True,
"font.family": "Helvetica"
})
def getPercentileLines(lineData,xs):
print(np.shape(lineData))
p05,p25,p50,p75,p95 = np.zeros(len(xs)), \
np.zeros(len(xs)), \
np.zeros(len(xs)), \
np.zeros(len(xs)), \
np.zeros(len(xs))
for i in range(len(xs)):
a = lineData[:,i]
p05[i] = np.percentile(a, 5)
p25[i] = np.percentile(a, 25)
p50[i] = np.percentile(a, 50)
p75[i] = np.percentile(a, 75)
p95[i] = np.percentile(a, 95)
return p05, p25, p50, p75, p95
orange = '#E69F00'
blue = '#0072B2'
green = '#009E73'
red = '#D55E00'
green2 = '#004D40'
blue2 = '#1E88E5'
orange2 = '#FFC107'
red2 = '#D81B60'
colSimp=orange
colGalExt=green
colGal = 'k'
agnHs = np.genfromtxt('agn_h_nF9.dat')
logfs = np.linspace(-9.,-6.,15)
p05,p25,p50,p75,p95 = getPercentileLines(agnHs,logfs)
p05,p25,p50,p75,p95 = getPercentileLines(agnHs,logfs)
print(np.shape(agnHs))
fs = [10.**f for f in logfs]
fs = np.atleast_1d(fs)
plt.fill_between(logfs,p05,p95,alpha=0.3,color=colSimp)
plt.fill_between(logfs,p25,p75,alpha=0.3,color=colSimp)
plt.plot(logfs,p50,color=colSimp,ls='--')
nFreqs=9
galLinesDir = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f{0}/f{0}_plotting_things/'.format(nFreqs)
hdata = np.genfromtxt('{}/h_lines.dat'.format(galLinesDir))
print(np.shape(hdata))
p05,p25,p50,p75,p95 = getPercentileLines(hdata,logfs)
fs = [10.**f for f in logfs]
fs = np.atleast_1d(fs)
plt.fill_between(logfs,p05,p95,alpha=0.3,color=colGalExt)
plt.fill_between(logfs,p25,p75,alpha=0.3,color=colGalExt)
plt.plot(logfs,p50,color=colGalExt,ls='--')
plt.ylabel(r'$h_{\rm c}$')
plt.xlabel(r'$\log_{10} f\,({\rm Hz})$')
plt.tight_layout()
plt.savefig('hlines.png')
plt.show()
plotting/gal_plotChains.py 0 → 100644
View file @ 47c5d4c9
import numpy as np
import matplotlib.pyplot as plt
pathToResults = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f9/'
run='ptmcmc_f9_2'
name = '_f9_2'
chains = np.genfromtxt(pathToResults+run+'/output/chain_1.txt')
nparameters = 18
fig, axes = plt.subplots(nparameters, figsize=(10, 10), sharex=True)
labels = np.genfromtxt('pars.txt',dtype=str)
print(labels)
#labels = ["m", "b", "log(f)"]
for i in range(nparameters):
ax = axes[i]
ax.plot(chains[:,i], "k", alpha=0.7)
ax.set_xlim(0, len(chains))
ax.set_ylabel(labels[i],rotation=0)
ax.yaxis.set_label_coords(-0.1, 0.5)
axes[-1].set_xlabel("step number");
plt.savefig('{0}{1}/chains{2}.png'.format(pathToResults,run,name))
plt.show()
plotting/nhistogram.py 0 → 100644
View file @ 47c5d4c9
from matplotlib import pylab
import numpy as np
from math import ceil
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'font.size': 16})
plt.rcParams.update({
"text.usetex": True,
"font.family": "Helvetica"
})
def fpairtest(delta):
return 1./(delta+1.)*(1.-0.25**(delta+1.))
def getneffGAL(nFreqBins):
if nFreqBins==9:
path = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f9/ptmcmc_f{}_1/output/'.format(nFreqBins)
elif nFreqBins==5:
path = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f5/ptmcmc_2/output/'.format(nFreqBins)
elif nFreqBins==3:
path = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f3/ptmcmc_f{}_1/output/'.format(nFreqBins)
chains = np.genfromtxt('{}/chain_1.txt'.format(path))
burnin = int(0.25*len(chains))
n1 = chains[burnin:,0]-chains[burnin:,2]+np.log10(chains[burnin:,5]/chains[burnin:,9]/fpairtest(chains[burnin:,8]))+chains[burnin:,10]*np.log10(0.4/0.7)+(chains[burnin:,10]-chains[burnin:,6])*(11.-chains[burnin:,2])+11.
priorChains = np.genfromtxt('/home/ADF/middlehr/repositories/PTAInference/runs/galaxyModel_ext/prior.dat')
priorn1 = priorChains[:,0]-priorChains[:,2]+np.log10(priorChains[:,5]/priorChains[:,9]/fpairtest(priorChains[:,8]))+priorChains[:,10]*np.log10(0.4/0.7)+(priorChains[:,10]-priorChains[:,6])*(11.-priorChains[:,2])+11.
return n1, priorn1
def getnAGN(nFreqBins):
if nFreqBins==9:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq9_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
elif nFreqBins==5:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq5_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
elif nFreqBins==3:
path0 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3/posterior.dat'.format(nFreqBins)
path1 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3_1/posterior.dat'.format(nFreqBins)
path2 = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{}/nFreq3_2/posterior.dat'.format(nFreqBins)
post0 = np.genfromtxt(path0,names=True)
post1 = np.genfromtxt(path1,names=True)
post2 = np.genfromtxt(path2,names=True)
#combine the ns
post0And1 = np.concatenate((post0['logn'], post1['logn']))
logns = np.concatenate((post0And1, post2['logn']))
return logns
def histOutline(dataIn, *args, **kwargs):
"""
code copied from http://www.scipy.org/Cookbook/Matplotlib/UnfilledHistograms?action=AttachFile&do=get&target=histNofill.py
Make a histogram that can be plotted with plot() so that
the histogram just has the outline rather than bars as it
usually does.
Example Usage:
binsIn = numpy.arange(0, 1, 0.1)
angle = pylab.rand(50)
(bins, data) = histOutline(binsIn, angle)
plot(bins, data, 'k-', linewidth=2)
"""
(histIn, binsIn) = np.histogram(dataIn, *args, **kwargs)
stepSize = binsIn[1] - binsIn[0]
bins = np.zeros(len(binsIn)*2 + 2, dtype=float)
data = np.zeros(len(binsIn)*2 + 2, dtype=float)
for bb in range(len(binsIn)):
bins[2*bb + 1] = binsIn[bb]
bins[2*bb + 2] = binsIn[bb] + stepSize
if bb < len(histIn):
data[2*bb + 1] = histIn[bb]
data[2*bb + 2] = histIn[bb]
bins[0] = bins[1]
bins[-1] = bins[-2]
data[0] = 0
data[-1] = 0
return (bins, data)
def plotOneDPosterior(data,paramLabel,setBins,colour='k-',label='',prior=False):
#pathToRuns = "../../runs/simpleModelPosteriors/"
#setBins = np.arange(int(min(data[paramName])), ceil(max(data[paramName])),0.5)
#setBins = np.linspace(min(data[paramName]), max(data[paramName]),50)
#setBins = np.linspace(-14.5,2.5,50)
p05 = np.percentile(data,5)
p50 = np.percentile(data,50)
p95 = np.percentile(data,95)
if prior==True:
plt.hist(data, bins=setBins, histtype='step', facecolor=colour, alpha=0.3, edgecolor=None,density=1,fill=True,label=label)
else:
(bins,dat) = histOutline(data,setBins,density=1)
pylab.plot(bins,dat,colour,linewidth=2,alpha=1.0,label=label)
plt.axvline(p05,ls=':',color=colour)
plt.axvline(p95,ls=':',color=colour)
plt.xlabel(paramLabel)
plt.ylabel(r'${\rm Normalised~counts}$')
plt.tight_layout()
#plt.savefig('{}lognoHist.pdf'.format(runLoc))
#plt.savefig('logno.png')
#pylab.show()
#print("""{}
#05%: {:.2f} ({})
#50%: {:.2f} ({})
#95%: {:.2f} ({})
#""".format(label,p05,10**p05,p50,10.**p50,p95,10.**p95))
#
print("""{}
{:.2f}^{{+{:.2f}}}_{{-{:.2f}}}""".format(label, p50, p95-p50, p50-p05))
return 0
n = 9
print('Number of frequency bins = {}'.format(n))
neffGAL, neffPriorGAL = getneffGAL(n)
nAGN = getnAGN(n)
# bins
low = int(min(nAGN))
high = ceil(max(nAGN))
bins = np.linspace(low,high,80)
orange = '#E69F00'
blue = '#0072B2'
green = '#009E73'
red = '#D55E00'
purple = '#CC79A7'
colSimp=orange
colGalExt=green
#colGal = 'k'
#colSimpNegAlpha = red
lognoLabel = r'$\log_{10}\frac{{\dot n}_0}{{\rm Mpc}^{3}{\rm Gyr}}$'
plotOneDPosterior(nAGN,lognoLabel,bins,colour=colSimp,label=r'$\textrm{Agnostic posterior}$')
plotOneDPosterior(neffGAL,lognoLabel,bins,colour=colGalExt,label=r'$\textrm{Astro-informed~posterior}$')
plotOneDPosterior(neffPriorGAL,lognoLabel,bins,colour=colGalExt,label=r'$\textrm{Astro-informed prior}$',prior=True)
plt.legend(fontsize=12)
plt.xlim(min(bins),max(bins))
plt.savefig('forPaper/nhistogram_nFreq{}.png'.format(n))
plt.savefig('forPaper/nhistogram_nFreq{}.pdf'.format(n))
plt.show()
plotting/pars.txt 0 → 100644
View file @ 47c5d4c9
Phi0
PhiI
M0
alpha0
alphaI
f0
betaf
alphaf
gammaf
tau0
alphatau
betatau
gammatau
Mstar
alphastar
epsilon
e0
zeta0
lnprob
lnlike
acc_rate
temp
plotting/plotData.py 0 → 100644
View file @ 47c5d4c9
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde as kde
def readChains(nFreqBins):
'''
read in the data and return the lowest nFreqBins
chains and frequencies
'''
frequencies = np.genfromtxt('../chains/dr2new_fs/freqs.txt')
samples = np.genfromtxt('../chains/dr2new_fs/chains_h.txt')
chainsToUse = samples[:,0:nFreqBins]
frequenciesToUse = frequencies[0:nFreqBins]
return frequenciesToUse, chainsToUse
def makeKDEs(samples,nFreqBins):
'''
make the KDEs to use in the likelihood
'''
KDE = {i: kde(samples[:,i]) for i in range(nFreqBins)}
return KDE
f,data = readChains(5)
kdes = makeKDEs(data,len(f))
hcsToPlot = np.linspace(-17,-13,20)
for i in range(len(f)):
chain = data[:,i]
kdeToPlot = kdes[i].logpdf(hcsToPlot)
print(kdeToPlot)
plt.hist(chain,density=True,histtype='step',bins=50)
plt.plot(hcsToPlot,kdeToPlot)
plt.show()
plotting/plotMLines.py 0 → 100644
View file @ 47c5d4c9
"""
plot distribution in chirp mass
used to make fugure 3 in paper
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'font.size': 16})
plt.rcParams.update({
"text.usetex": True,
"font.family": "Helvetica"
})
def getPercentileLines(lineData,xs):
print(np.shape(lineData))
p01,p05,p25,p50,p75,p95,p99 = np.zeros(len(xs)), \
np.zeros(len(xs)), \
np.zeros(len(xs)), \
np.zeros(len(xs)), \
np.zeros(len(xs)), \
np.zeros(len(xs)), \
np.zeros(len(xs))
for i in range(len(xs)):
a = lineData[:,i]
p01[i] = np.percentile(a, 0.5)
p05[i] = np.percentile(a, 5)
p25[i] = np.percentile(a, 25)
p50[i] = np.percentile(a, 50)
p75[i] = np.percentile(a, 75)
p95[i] = np.percentile(a, 95)
p99[i] = np.percentile(a, 99.5)
return p01, p05, p25, p50, p75, p95, p99
orange = '#E69F00'
blue = '#0072B2'
green = '#009E73'
red = '#D55E00'
green2 = '#004D40'
blue2 = '#1E88E5'
orange2 = '#FFC107'
red2 = '#D81B60'
colSimp=orange
colGalExt=green
colGal = 'k'
nFreqs = 9
agnLinesDir = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/agn_model/f{0}/f{0}_plotting_things/'.format(nFreqs)
lines = np.genfromtxt('{}/mLines.dat'.format(agnLinesDir))
simpModMs = np.genfromtxt('{}/ms.dat'.format(agnLinesDir))
print(type(simpModMs),np.shape(simpModMs))
_, p05, p25, p50, p75, p95, _ = getPercentileLines(lines,simpModMs)
plt.fill_between(simpModMs,p05,p95,alpha=0.3,color=colSimp)
plt.fill_between(simpModMs,p25,p75,alpha=0.3,color=colSimp)
plt.plot(simpModMs,p50,color=colSimp,ls='--')
plt.ylim(1E-9,1E6)
plt.xlim(1E6,1E11)
plt.yscale('log')
plt.xscale('log')
plt.xlabel(r'$\mathcal{M}~~({\rm M_{\odot}})$')
plt.ylabel(r'${\rm d}n / {\rm d} \log_{10} \mathcal{M}/{\rm M_{\odot}}~~({\rm Mpc}^{-3})$')
plt.tight_layout()
#plt.savefig('dndlogM_3FreqBins.png')
galLinesDir = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f{0}/f{0}_plotting_things/'.format(nFreqs)
galExtModData = np.genfromtxt('{}/m_lines.dat'.format(galLinesDir))
galExtModLogMs = np.linspace(5,11,30)
print(type(galExtModData),np.shape(galExtModData))
galExtModMs = [ 10.0**logM for logM in galExtModLogMs ]
galExtModMs = np.atleast_1d(galExtModMs)
_, p05,p25,p50,p75,p95, _ = getPercentileLines(galExtModData,galExtModMs)
plt.fill_between(galExtModMs,p05,p95,alpha=0.3,color=colGalExt)
plt.fill_between(galExtModMs,p25,p75,alpha=0.3,color=colGalExt,label='new posterior')
plt.plot(galExtModMs,p50,color=colGalExt,ls='--')
#plt.savefig('forPaper/mlines_nFreq{}.pdf'.format(nFreqs))
#plt.savefig('forPaper/mlines_nFreq{}.png'.format(nFreqs))
#plt.show()
#exit()
galPriorLinesDir = '/home/ADF/middlehr/Documents/PTA/EPTAInterp/results/gal_model/f{0}/f{0}_plotting_things_prior/'.format(nFreqs)
galExtModPrior = np.genfromtxt('{}/m_linesprior.dat'.format(galPriorLinesDir))
p01,p05,p25,p50,p75,p95,p99= getPercentileLines(galExtModPrior,galExtModMs)
plt.plot(galExtModMs, p01, color='k', alpha=0.7, ls=':',label='newly computed prior')
plt.plot(galExtModMs, p99, color='k', alpha=0.7, ls=':')
"""
galModRunLoc = '../../runs/galaxyModel/'
galModData = np.genfromtxt('{}/m_lines.dat'.format(galModRunLoc))
galModLogMs = np.genfromtxt('{}/mc.txt'.format(galModRunLoc))
galModMs = [ 10.0**logM for logM in galModLogMs ]
p05,p25,p50,p75,p95 = confBands.getPercentiles(galModMs,galModData)
plt.fill_between(galModMs,p05,p95,alpha=0.3,color=colGal)
plt.fill_between(galModMs,p25,p75,alpha=0.3,color=colGal)
plt.plot(galModMs,p50,color=colGal,ls='--')
"""
galPriorLinesDir = '/home/ADF/middlehr/repositories/PTAInference/runs/galaxyModel_ext/'
galExtModPrior = np.genfromtxt('{}/m_linesprior.dat'.format(galPriorLinesDir))
p01,p05,p25,p50,p75,p95,p99= getPercentileLines(galExtModPrior,galExtModMs)
plt.plot(galExtModMs, p01, color='k', alpha=0.7, ls='-.',label='M21 prior')
plt.plot(galExtModMs, p99, color='k', alpha=0.7, ls='-.')
galLinesDir = '/home/ADF/middlehr/repositories/PTAInference/runs/galaxyModel_ext/'
galExtModData = np.genfromtxt('{}/m_lines.dat'.format(galLinesDir))
_, p05,p25,p50,p75,p95, _ = getPercentileLines(galExtModData,galExtModMs)
plt.fill_between(galExtModMs,p05,p95,alpha=0.2,color='b',label='M21 posterior')
plt.fill_between(galExtModMs,p25,p75,alpha=0.2,color='b')
plt.plot(galExtModMs,p50,color='b',ls='--')
plt.legend()
plt.ylim(1E-9,1E6)
plt.xlim(1E6,1E11)
plt.yscale('log')
plt.xscale('log')
plt.xlabel(r'$\mathcal{M}~~({\rm M_{\odot}})$')
plt.ylabel(r'${\rm d}n / {\rm d} \log_{10} \mathcal{M}/{\rm M_{\odot}}~~({\rm Mpc}^{-3})$')
plt.tight_layout()
plt.savefig('forPaper/mlines_nFreq{}.pdf'.format(nFreqs))
plt.savefig('forPaper/mlines_nFreq{}.png'.format(nFreqs))
plt.show()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment