Commit d3b8d4d8 authored by Brandon Piotrzkowski's avatar Brandon Piotrzkowski

Create unit tests for priors.py

parent 1b51472e
Pipeline #83704 passed with stages
in 2 minutes and 21 seconds
from gwcosmo.prior import priors
import numpy as np
from scipy import integrate
import math as mth
def test_pH0():
assert priors.pH0(70) == 1/70
assert priors.pH0(70, prior='log') == 1/70
assert priors.pH0([70], prior='uniform') == 1.
def test_BBH_mass_distribution():
m_min = 4.
m_max = 40.
alpha = 1.6
n = priors.BBH_mass_distribution(1000000,
mmin=m_min,
mmax=m_max,
alpha=alpha)
# Define analytic prob density funcitons
def p(m, alpha):
return m**(-alpha)
def norm(m_min, m_max, alpha):
return integrate.quad(p, m_min, m_max, args=(alpha))[0]
def prob(m, mmin, mmax, alpha):
return p(m, alpha)/norm(m_min, m_max, alpha)
# Bin up counts from above and get prob density
counts, masses = np.histogram(n[0], bins=100, density=True)
# Get masses from center of bin
masses = (masses[0:-1] + masses[1:])/2
# Get analytic results to compare to
probs = []
for mass in masses:
probs.append(prob(mass, m_min, m_max, alpha))
# Compare numerical and analytic results
is_close = []
for i in np.arange(len(counts)):
is_close.append(mth.isclose(counts[i],probs[i], rel_tol=.05))
# Assert that a number of bins should be within tolerance
print(sum(is_close)/len(is_close))
assert sum(is_close)/len(is_close) > .80
def test_BNS_guassian_distribution():
m_min = .6
m_max = 2.1
mean = 1.35
sigma = 0.15
m1, m2 = priors.BNS_gaussian_distribution(1000000,
mean=mean,
sigma=sigma)
# Define analytic prob density funcitons
def p(m, mean, sigma):
return 1/(2*np.pi*sigma**2)**(.5)*np.exp(-(m-mean)**2./(2.*sigma**2.))
def norm(m_min, m_max, mean, sigma):
return integrate.quad(p, m_min, m_max, args=(mean, sigma))[0]
def prob(m, mmin, mmax, mean, sigma):
return p(m, mean, sigma)/norm(m_min, m_max, mean, sigma)
# Bin up counts from above and get prob density
counts, masses = np.histogram(np.concatenate((m1,m2)), bins=100, density=True)
# Get masses from center of bin
masses = (masses[0:-1] + masses[1:])/2
# Get analytic results to compare to
probs = []
for mass in masses:
probs.append(prob(mass, m_min, m_max, mean, sigma))
# Compare numerical and analytic results
is_close = []
for i in np.arange(len(counts)):
is_close.append(mth.isclose(counts[i],probs[i], abs_tol=.01))
# Assert that a number of bins should be within tolerance
print(sum(is_close)/len(is_close))
assert sum(is_close)/len(is_close) > .80
def test_BNS_uniform_distribution():
m_min = 1.
m_max = 3.
m1, m2 = priors.BNS_uniform_distribution(1000000,
mmin=m_min,
mmax=m_max)
# Define analytic prob density funcitons
def prob(mmin, mmax):
return 1/(mmax-mmin)
# Bin up counts from above and get prob density
counts, masses = np.histogram(np.concatenate((m1,m2)), bins=100, density=True)
# Get masses from center of bin
masses = (masses[0:-1] + masses[1:])/2
# Get analytic results to compare to
probs = []
for mass in masses:
probs.append(prob(m_min, m_max))
# Compare numerical and analytic results
is_close = []
for i in np.arange(len(counts)):
is_close.append(mth.isclose(counts[i],probs[i], abs_tol=.01))
# Assert that a number of bins should be within tolerance
print(sum(is_close)/len(is_close))
assert sum(is_close)/len(is_close) > .80
......@@ -13,3 +13,4 @@ matplotlib>=3.1.1
numpy>=1.16.4
pandas>=0.24.2
progressbar
scipy
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