Merge branch 'master' of http://gitlab.sr.bham.ac.uk/hannahm/PTAInterpretation

models/galaxy_model/mergerrate.py

@@ -23,12 +23,12 @@ def dtdz(z):
 
 def dVcdz(z):
 	c = const.c.value  # 3.e8
-	H0 = cosmo.H(0).value*1.e3 # 70.e3 -> check with Siyuan
+	H0 = cosmo.H(0).value*1.e3 # 70.e3 
 	return 4.*np.pi*c/H0*invE(z)*DM(z)*DM(z)
 
 def DM(z):
 	c = const.c.value   # 3.e8
-	H0 = cosmo.H(0).value*1.e3 # 70.e3 -> check with Siyuan
+	H0 = cosmo.H(0).value*1.e3 # 70.e3 
 	return c/H0*si.quad(invE, 0., z)[0]
 
 def mchirpq(q):
@@ -53,7 +53,7 @@ def mfractionp(m):
 #spectrum computation functions begin
 
 Gu = const.G.value                # 6.673848e-11
-cu = const.G.value                # 2.997925e8
+cu = const.c.value                # 2.997925e8
 Msun = const.M_sun.value          # 1.98855e30
 pc = 1.*u.parsec.to(u.m)          # 3.0856776e16
 Gyr = 1.e9 * 1.*u.year.to(u.s)    # 3.15576e16
@@ -404,8 +404,8 @@ class mergerrate(object):
 	def hdrop(self,n0,f,fbin):
 		if fbin is None:
 			return n0
-		c = const.c.value   # 3.e8
-		G = 1.33e20 # should we use const.G.value here too (not sure what the unit is here)
+		c = cu   # 3.e8
+		G = const.G.value*const.M_sun.value  # 1.33e20 
 		fhigh = f + fbin
 		flow = f - fbin
 		Mcbh = np.linspace(5,11,30)

models/galaxy_model/run_sampler.py

+# This is based on the ptmcmcrun.py here:
+# https://github.com/hannahm8/PTAInference/tree/master/code/galaxy_model 
+# 
+
+import numpy as np
+from scipy.stats import gaussian_kde as kde
+import matplotlib.pyplot as plt
+
+
+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(nFrequencies)}
+
+    return KDE
+
+
+def in_bounds(par):
+    return all(bounds[i][0] < par[i] < bounds[i][1] for i in range(len(par)))
+
+
+def log_prior(par):
+    if in_bounds(par):
+        param = dict(Phi0 = par[0], PhiI = par[1], M0 = par[2], alpha = par[3], alphaI = par[4], f0 = par[5], beta = par[6], gamma = par[7], delta = par[8], t0 = par[9], epsilon = par[10], zeta = par[11], eta = par[12], Ms = par[13], theta = par[14], sigma = par[15], e0 = par[16], rho = par[17])        
+        return prior_check.check_p(param)
+    else:
+        return -np.inf
+
+
+def log_likelihood(par):
+    
+    M1 = np.linspace(9,12,25)
+    q = np.linspace(0.25,1,10)
+    z = np.linspace(0.,1.5,5)
+
+    initpar = dict(Phi0 = par[0], 
+                   PhiI = par[1], 
+                   M0 = par[2], 
+                   alpha = par[3], 
+                   alphaI = par[4], 
+                   f0 = par[5], 
+                   beta = par[6], 
+                   gamma = par[7], 
+                   delta = par[8], 
+                   t0 = par[9], 
+                   epsilon = par[10], 
+                   zeta = par[11], 
+                   eta = par[12], 
+                   Ms = par[13], 
+                   theta = par[14], 
+                   sigma = par[15], 
+                   e0 = par[16], 
+                   rho = par[17])
+    model = mr.mergerrate(M1,q,z,f,**initpar).hmodelt(fbin=None)[0]
+
+    logL = 0.
+    for i in range(nFrequencies):
+        logL += KDE[i].logpdf(model)
+    return logL
+
+
+
+
+numberOfFrequencyBinsToUse = 5
+frequencies, samples = readChains(numberOfFrequencyBinsToUse)
+KDEs = makeKDEs(samples, numberOfFrequencyBinsToUse)
+
+
+
+to do - add sampler
+
+
+
+'''
+some testing...
+KDE = {i: kde(samples[:,i]) for i in range(nFrequencies)}
+print(KDE)
+
+
+for i in range(nFrequencies):
+
+    bins = np.linspace(min(samples[:,i]), max(samples[:,i]), 100)
+    print('bins' , bins)
+    toPlot = np.linspace(min(samples[:,i]), max(samples[:,i]), 100)
+    print(KDE[i].logpdf(-14))
+    
+    plt.plot(toPlot,KDE[i].logpdf(toPlot))
+    plt.hist(samples[:,i],bins=bins,density=True,histtype='step')
+    plt.show()
+
+    
+'''