Skip to content
Snippets Groups Projects
Commit 0f0d24ff authored by Moritz Huebner's avatar Moritz Huebner
Browse files

Merge branch 'update_prior_files' into 'master'

Add units to prior files

See merge request Monash/tupak!214
parents 0019ed76 28d6dba7
No related merge requests found
Showing
with 108 additions and 76 deletions
......@@ -25,6 +25,7 @@ Changes currently on master, but not under a tag.
- Added NestedSampler and MCSampler helper classes
- Added sampler_requirements.txt file
- Add AlignedSpin gw prior
- Add units to know prior files
### Changes
- Fix construct_cbc_derived_parameters
......
......@@ -60,8 +60,9 @@ ifos.inject_signal(waveform_generator=waveform_generator,
# that will be included in the sampler. If we do nothing, then the default priors get used.
priors = tupak.gw.prior.BBHPriorSet()
priors['geocent_time'] = tupak.core.prior.Uniform(
minimum=injection_parameters['geocent_time'] - 1, maximum=injection_parameters['geocent_time'] + 1,
name='geocent_time', latex_label='$t_c$')
minimum=injection_parameters['geocent_time'] - 1,
maximum=injection_parameters['geocent_time'] + 1,
name='geocent_time', latex_label='$t_c$', unit='$s$')
for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi', 'ra', 'dec', 'geocent_time', 'phase']:
priors[key] = injection_parameters[key]
......
......@@ -50,7 +50,8 @@ priors.pop('mass_1')
priors.pop('mass_2')
priors.pop('luminosity_distance')
priors['chirp_mass'] = tupak.prior.Uniform(
name='chirp_mass', latex_label='$m_c$', minimum=13, maximum=45)
name='chirp_mass', latex_label='$m_c$', minimum=13, maximum=45,
unit='$M_{\\odot}$')
priors['symmetric_mass_ratio'] = tupak.prior.Uniform(
name='symmetric_mass_ratio', latex_label='q', minimum=0.1, maximum=0.25)
priors['redshift'] = tupak.prior.Uniform(
......
......@@ -47,8 +47,9 @@ ifos.inject_signal(waveform_generator=waveform,
# create the priors
prior = injection_parameters.copy()
prior['amplitude'] = tupak.core.prior.Uniform(1e-23, 1e-21, r'$h_0$')
prior['damping_time'] = tupak.core.prior.Uniform(0, 1, r'damping time')
prior['frequency'] = tupak.core.prior.Uniform(0, 200, r'frequency')
prior['damping_time'] = tupak.core.prior.Uniform(
0, 1, r'damping time', unit='$s$')
prior['frequency'] = tupak.core.prior.Uniform(0, 200, r'frequency', unit='Hz')
prior['phase'] = tupak.core.prior.Uniform(-np.pi / 2, np.pi / 2, r'$\phi$')
# define likelihood
......
......@@ -60,16 +60,23 @@ ifos.inject_signal(waveform_generator=waveform_generator,
# Now we set up the priors on each of the binary parameters.
priors = dict()
priors["mass_1"] = tupak.core.prior.Uniform(name='mass_1', minimum=5, maximum=60)
priors["mass_2"] = tupak.core.prior.Uniform(name='mass_2', minimum=5, maximum=60)
priors["eccentricity"] = tupak.core.prior.PowerLaw(name='eccentricity', latex_label='$e$', alpha=-1, minimum=1e-4, maximum=0.4)
priors["luminosity_distance"] = tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=2e3)
priors["dec"] = tupak.core.prior.Cosine(name='dec')
priors["ra"] = tupak.core.prior.Uniform(name='ra', minimum=0, maximum=2 * np.pi)
priors["iota"] = tupak.core.prior.Sine(name='iota')
priors["psi"] = tupak.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi)
priors["phase"] = tupak.core.prior.Uniform(name='phase', minimum=0, maximum=2 * np.pi)
priors["geocent_time"] = tupak.core.prior.Uniform(1180002600.9, 1180002601.1, name='geocent_time')
priors["mass_1"] = tupak.core.prior.Uniform(
name='mass_1', minimum=5, maximum=60, unit='$M_{\\odot}$')
priors["mass_2"] = tupak.core.prior.Uniform(
name='mass_2', minimum=5, maximum=60, unit='$M_{\\odot}$')
priors["eccentricity"] = tupak.core.prior.PowerLaw(
name='eccentricity', latex_label='$e$', alpha=-1, minimum=1e-4, maximum=0.4)
priors["luminosity_distance"] = tupak.gw.prior.UniformComovingVolume(
name='luminosity_distance', minimum=1e2, maximum=2e3)
priors["dec"] = tupak.core.prior.Cosine(name='dec')
priors["ra"] = tupak.core.prior.Uniform(
name='ra', minimum=0, maximum=2 * np.pi)
priors["iota"] = tupak.core.prior.Sine(name='iota')
priors["psi"] = tupak.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi)
priors["phase"] = tupak.core.prior.Uniform(
name='phase', minimum=0, maximum=2 * np.pi)
priors["geocent_time"] = tupak.core.prior.Uniform(
1180002600.9, 1180002601.1, name='geocent_time', unit='s')
# Initialising the likelihood function.
likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers=ifos,
......
......@@ -40,25 +40,29 @@ ifos.inject_signal(waveform_generator=waveform_generator,
# This loads in a predefined set of priors for BBHs.
priors = tupak.gw.prior.BBHPriorSet()
# These parameters will not be sampled
for key in ['tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'iota', 'ra', 'dec', 'geocent_time', 'psi']:
for key in ['tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'iota', 'ra',
'dec', 'geocent_time', 'psi']:
priors[key] = injection_parameters[key]
# We can make uniform distributions.
priors['mass_2'] = tupak.core.prior.Uniform(name='mass_2', minimum=20, maximum=40)
priors['mass_2'] = tupak.core.prior.Uniform(
name='mass_2', minimum=20, maximum=40, unit='$M_{\\odot}$')
# We can make a power-law distribution, p(x) ~ x^{alpha}
# Note: alpha=0 is a uniform distribution, alpha=-1 is uniform-in-log
priors['a_1'] = tupak.core.prior.PowerLaw(name='a_1', alpha=-1, minimum=1e-2, maximum=1)
priors['a_1'] = tupak.core.prior.PowerLaw(
name='a_1', alpha=-1, minimum=1e-2, maximum=1)
# We can define a prior from an array as follows.
# Note: this doesn't have to be properly normalised.
a_2 = np.linspace(0, 1, 1001)
p_a_2 = a_2 ** 4
priors['a_2'] = tupak.core.prior.Interped(name='a_2', xx=a_2, yy=p_a_2, minimum=0, maximum=0.5)
priors['a_2'] = tupak.core.prior.Interped(
name='a_2', xx=a_2, yy=p_a_2, minimum=0, maximum=0.5)
# Additionally, we have Gaussian, TruncatedGaussian, Sine and Cosine.
# It's also possible to load an interpolate a prior from a file.
# Finally, if you don't specify any necessary parameters it will be filled in from the default when the sampler starts.
# Enjoy.
# Initialise GravitationalWaveTransient
likelihood = tupak.GravitationalWaveTransient(interferometers=ifos, waveform_generator=waveform_generator)
likelihood = tupak.gw.GravitationalWaveTransient(interferometers=ifos, waveform_generator=waveform_generator)
# Run sampler
result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty',
......
......@@ -53,7 +53,7 @@ for key in ['psi', 'ra', 'dec', 'geocent_time']:
#priors['Q'] = tupak.prior.create_default_prior(name='Q')
#priors['frequency'] = tupak.prior.create_default_prior(name='frequency')
priors['Q'] = tupak.core.prior.Uniform(2, 50, 'Q')
priors['frequency'] = tupak.core.prior.Uniform(30, 1000, 'frequency')
priors['frequency'] = tupak.core.prior.Uniform(30, 1000, 'frequency', unit='Hz')
priors['hrss'] = tupak.core.prior.Uniform(1e-23, 1e-21, 'hrss')
# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
......
......@@ -68,7 +68,7 @@ priors = dict()
for key in ['psi', 'geocent_time']:
priors[key] = injection_parameters[key]
priors['luminosity_distance'] = tupak.core.prior.Uniform(
2, 20, 'luminosity_distance')
2, 20, 'luminosity_distance', unit='$kpc$')
priors['pc_coeff1'] = tupak.core.prior.Uniform(-1, 1, 'pc_coeff1')
priors['pc_coeff2'] = tupak.core.prior.Uniform(-1, 1, 'pc_coeff2')
priors['pc_coeff3'] = tupak.core.prior.Uniform(-1, 1, 'pc_coeff3')
......@@ -80,7 +80,7 @@ priors['dec'] = tupak.core.prior.Sine(name='dec')
priors['geocent_time'] = tupak.core.prior.Uniform(
injection_parameters['geocent_time'] - 1,
injection_parameters['geocent_time'] + 1,
'geocent_time')
'geocent_time', unit='$s$')
# Initialise the likelihood by passing in the interferometer data (IFOs) and
# the waveoform generator
......
......@@ -2,10 +2,10 @@
# Note that you may wish to use more specific mass and distance parameters.
# These commands are all known to tupak.gw.prior.
# Lines beginning "#" are ignored.
mass_1 = Uniform(name='mass_1', minimum=5, maximum=100)
mass_2 = Uniform(name='mass_2', minimum=5, maximum=100)
# chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=100)
# total_mass = Uniform(name='total_mass', minimum=10, maximum=200)
mass_1 = Uniform(name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$')
mass_2 = Uniform(name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$')
# chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=100, unit='$M_{\\odot}$')
# total_mass = Uniform(name='total_mass', minimum=10, maximum=200, unit='$M_{\\odot}$')
# mass_ratio = Uniform(name='mass_ratio', minimum=0.125, maximum=1)
# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=8 / 81, maximum=0.25)
a_1 = Uniform(name='a_1', minimum=0, maximum=0.8)
......@@ -16,7 +16,7 @@ tilt_2 = Sine(name='tilt_2')
# cos_tilt_2 = Uniform(name='cos_tilt_2', minimum=-1, maximum=1)
phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi)
phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
luminosity_distance = tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3)
luminosity_distance = tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc')
dec = Cosine(name='dec')
ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
iota = Sine(name='iota')
......
......@@ -347,22 +347,30 @@ class TestPriorSet(unittest.TestCase):
# shutil.rmtree(outdir_path)
def test_read_from_file(self):
expected = dict(mass_1=tupak.core.prior.Uniform(name='mass_1', minimum=5, maximum=100),
mass_2=tupak.core.prior.Uniform(name='mass_2', minimum=5, maximum=100),
a_1=tupak.core.prior.Uniform(name='a_1', minimum=0, maximum=0.8),
a_2=tupak.core.prior.Uniform(name='a_2', minimum=0, maximum=0.8),
tilt_1=tupak.core.prior.Sine(name='tilt_1'),
tilt_2=tupak.core.prior.Sine(name='tilt_2'),
phi_12=tupak.core.prior.Uniform(name='phi_12', minimum=0, maximum=2 * np.pi),
phi_jl=tupak.core.prior.Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi),
luminosity_distance=tupak.gw.prior.UniformComovingVolume(name='luminosity_distance',
minimum=1e2, maximum=5e3),
dec=tupak.core.prior.Cosine(name='dec'),
ra=tupak.core.prior.Uniform(name='ra', minimum=0, maximum=2 * np.pi),
iota=tupak.core.prior.Sine(name='iota'),
psi=tupak.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi),
phase=tupak.core.prior.Uniform(name='phase', minimum=0, maximum=2 * np.pi)
)
expected = dict(
mass_1=tupak.core.prior.Uniform(
name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$'),
mass_2=tupak.core.prior.Uniform(
name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$'),
a_1=tupak.core.prior.Uniform(name='a_1', minimum=0, maximum=0.8),
a_2=tupak.core.prior.Uniform(name='a_2', minimum=0, maximum=0.8),
tilt_1=tupak.core.prior.Sine(name='tilt_1'),
tilt_2=tupak.core.prior.Sine(name='tilt_2'),
phi_12=tupak.core.prior.Uniform(
name='phi_12', minimum=0, maximum=2 * np.pi),
phi_jl=tupak.core.prior.Uniform(
name='phi_jl', minimum=0, maximum=2 * np.pi),
luminosity_distance=tupak.gw.prior.UniformComovingVolume(
name='luminosity_distance', minimum=1e2,
maximum=5e3, unit='Mpc'),
dec=tupak.core.prior.Cosine(name='dec'),
ra=tupak.core.prior.Uniform(
name='ra', minimum=0, maximum=2 * np.pi),
iota=tupak.core.prior.Sine(name='iota'),
psi=tupak.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi),
phase=tupak.core.prior.Uniform(
name='phase', minimum=0, maximum=2 * np.pi)
)
self.assertDictEqual(expected, self.prior_set_from_file)
def test_convert_floats_to_delta_functions(self):
......@@ -380,22 +388,30 @@ class TestPriorSet(unittest.TestCase):
def test_prior_set_from_dict_but_using_a_string(self):
prior_set = tupak.core.prior.PriorSet(dictionary=self.default_prior_file)
expected = dict(mass_1=tupak.core.prior.Uniform(name='mass_1', minimum=5, maximum=100),
mass_2=tupak.core.prior.Uniform(name='mass_2', minimum=5, maximum=100),
a_1=tupak.core.prior.Uniform(name='a_1', minimum=0, maximum=0.8),
a_2=tupak.core.prior.Uniform(name='a_2', minimum=0, maximum=0.8),
tilt_1=tupak.core.prior.Sine(name='tilt_1'),
tilt_2=tupak.core.prior.Sine(name='tilt_2'),
phi_12=tupak.core.prior.Uniform(name='phi_12', minimum=0, maximum=2 * np.pi),
phi_jl=tupak.core.prior.Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi),
luminosity_distance=tupak.gw.prior.UniformComovingVolume(name='luminosity_distance',
minimum=1e2, maximum=5e3),
dec=tupak.core.prior.Cosine(name='dec'),
ra=tupak.core.prior.Uniform(name='ra', minimum=0, maximum=2 * np.pi),
iota=tupak.core.prior.Sine(name='iota'),
psi=tupak.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi),
phase=tupak.core.prior.Uniform(name='phase', minimum=0, maximum=2 * np.pi)
)
expected = dict(
mass_1=tupak.core.prior.Uniform(
name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$'),
mass_2=tupak.core.prior.Uniform(
name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$'),
a_1=tupak.core.prior.Uniform(name='a_1', minimum=0, maximum=0.8),
a_2=tupak.core.prior.Uniform(name='a_2', minimum=0, maximum=0.8),
tilt_1=tupak.core.prior.Sine(name='tilt_1'),
tilt_2=tupak.core.prior.Sine(name='tilt_2'),
phi_12=tupak.core.prior.Uniform(
name='phi_12', minimum=0, maximum=2 * np.pi),
phi_jl=tupak.core.prior.Uniform(
name='phi_jl', minimum=0, maximum=2 * np.pi),
luminosity_distance=tupak.gw.prior.UniformComovingVolume(
name='luminosity_distance', minimum=1e2,
maximum=5e3, unit='Mpc'),
dec=tupak.core.prior.Cosine(name='dec'),
ra=tupak.core.prior.Uniform(
name='ra', minimum=0, maximum=2 * np.pi),
iota=tupak.core.prior.Sine(name='iota'),
psi=tupak.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi),
phase=tupak.core.prior.Uniform(
name='phase', minimum=0, maximum=2 * np.pi)
)
self.assertDictEqual(expected, prior_set)
def test_dict_argument_is_not_string_or_dict(self):
......
......@@ -8,7 +8,8 @@ from scipy.interpolate import UnivariateSpline
class UniformComovingVolume(FromFile):
def __init__(self, minimum=None, maximum=None, name='luminosity distance', latex_label='$d_L$'):
def __init__(self, minimum=None, maximum=None,
name='luminosity distance', latex_label='$d_L$', unit='Mpc'):
"""
Parameters
......@@ -24,7 +25,7 @@ class UniformComovingVolume(FromFile):
"""
file_name = os.path.join(os.path.dirname(__file__), 'prior_files', 'comoving.txt')
FromFile.__init__(self, file_name=file_name, minimum=minimum, maximum=maximum, name=name,
latex_label=latex_label, unit='Mpc')
latex_label=latex_label, unit=unit)
class AlignedSpin(Interped):
......
# These are the default priors for analysing GW150914.
mass_1 = Uniform(name='mass_1', minimum=30, maximum=50)
mass_2 = Uniform(name='mass_2', minimum=20, maximum=40)
mass_1 = Uniform(name='mass_1', minimum=30, maximum=50, unit='$M_{\\odot}$')
mass_2 = Uniform(name='mass_2', minimum=20, maximum=40, unit='$M_{\\odot}$')
a_1 = Uniform(name='a_1', minimum=0, maximum=0.8)
a_2 = Uniform(name='a_2', minimum=0, maximum=0.8)
tilt_1 = Sine(name='tilt_1')
tilt_2 = Sine(name='tilt_2')
phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi)
phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
luminosity_distance = tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=1e3)
luminosity_distance = tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=1e3, unit='Mpc')
dec = Cosine(name='dec')
ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
iota = Sine(name='iota')
psi = Uniform(name='psi', minimum=0, maximum=np.pi)
phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi)
geocent_time = Uniform(1126259462.322, 1126259462.522, name='geocent_time')
geocent_time = Uniform(1126259462.322, 1126259462.522, name='geocent_time', unit='$s$')
# These are the calibration parameters as described in
# https://journals.aps.org/prx/abstract/10.1103/PhysRevX.6.041015
# recalib_H1_frequency_0 = 20
......
......@@ -2,10 +2,10 @@
# Note that you may wish to use more specific mass and distance parameters.
# These commands are all known to tupak.gw.prior.
# Lines beginning "#" are ignored.
mass_1 = Uniform(name='mass_1', minimum=5, maximum=100)
mass_2 = Uniform(name='mass_2', minimum=5, maximum=100)
# chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=100)
# total_mass = Uniform(name='total_mass', minimum=10, maximum=200)
mass_1 = Uniform(name='mass_1', minimum=5, maximum=100, unit='$M_{\\odot}$')
mass_2 = Uniform(name='mass_2', minimum=5, maximum=100, unit='$M_{\\odot}$')
# chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=100, unit='$M_{\\odot}$')
# total_mass = Uniform(name='total_mass', minimum=10, maximum=200, unit='$M_{\\odot}$')
# mass_ratio = Uniform(name='mass_ratio', minimum=0.125, maximum=1)
# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=8 / 81, maximum=0.25)
a_1 = Uniform(name='a_1', minimum=0, maximum=0.8)
......@@ -16,7 +16,7 @@ tilt_2 = Sine(name='tilt_2')
# cos_tilt_2 = Uniform(name='cos_tilt_2', minimum=-1, maximum=1)
phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi)
phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
luminosity_distance = tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3)
luminosity_distance = tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc')
dec = Cosine(name='dec')
ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
iota = Sine(name='iota')
......
......@@ -2,15 +2,15 @@
# Note that you may wish to use more specific mass and distance parameters.
# These commands are all known to tupak.gw.prior.
# Lines beginning "#" are ignored.
mass_1 = Uniform(name='mass_1', minimum=1, maximum=2)
mass_2 = Uniform(name='mass_2', minimum=1, maximum=2)
# chirp_mass = Uniform(name='chirp_mass', minimum=0.87, maximum=1.74)
# total_mass = Uniform(name='total_mass', minimum=2, maximum=4)
mass_1 = Uniform(name='mass_1', minimum=1, maximum=2, unit='$M_{\\odot}$')
mass_2 = Uniform(name='mass_2', minimum=1, maximum=2, unit='$M_{\\odot}$')
# chirp_mass = Uniform(name='chirp_mass', minimum=0.87, maximum=1.74, unit='$M_{\\odot}$')
# total_mass = Uniform(name='total_mass', minimum=2, maximum=4, unit='$M_{\\odot}$')
# mass_ratio = Uniform(name='mass_ratio', minimum=0.5, maximum=1)
# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=0.22, maximum=0.25)
a_1 = Uniform(name='a_1', minimum= -0.05, maximum= 0.05)
a_2 = Uniform(name='a_2', minimum= -0.05, maximum= 0.05)
luminosity_distance = tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500)
luminosity_distance = tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500, unit='Mpc')
dec = Cosine(name='dec')
ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
iota = Sine(name='iota')
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment