From ff805de92376bea5dc42dfd74b47fb0df4e461f2 Mon Sep 17 00:00:00 2001
From: Gregory Ashton <gregory.ashton@ligo.org>
Date: Sun, 10 Jun 2018 16:35:14 +1000
Subject: [PATCH] Major fix up to priors
- Moves create_default_prior to the core prior file
- Makes fill_prior a class attribute
- Makes the BBH default priors its own file which is read in when
figuring out the defaults
- Checks the example all still run
---
examples/injection_examples/basic_tutorial.py | 2 +-
.../basic_tutorial_dist_time_phase_marg.py | 2 +-
.../basic_tutorial_time_phase_marg.py | 2 +-
.../change_sampled_parameters.py | 2 +-
.../how_to_specify_the_prior.py | 2 +-
.../supernova_example/supernova_example.py | 4 +-
tupak/core/prior.py | 149 +++++++++++-------
.../core/prior_files/binary_black_holes.prior | 21 +++
tupak/core/sampler.py | 4 +-
tupak/gw/likelihood.py | 4 +-
tupak/gw/prior.py | 52 +-----
11 files changed, 124 insertions(+), 120 deletions(-)
create mode 100644 tupak/core/prior_files/binary_black_holes.prior
diff --git a/examples/injection_examples/basic_tutorial.py b/examples/injection_examples/basic_tutorial.py
index 85cb46977..00e8403b2 100644
--- a/examples/injection_examples/basic_tutorial.py
+++ b/examples/injection_examples/basic_tutorial.py
@@ -56,7 +56,7 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl','psi', 'ra', 'd
# The above list does *not* include mass_1, mass_2, iota and luminosity_distance, which means those are the parameters
# that will be included in the sampler. If we do nothing, then the default priors get used.
-priors['luminosity_distance'] = tupak.gw.prior.create_default_prior(name='luminosity_distance')
+priors['luminosity_distance'] = tupak.core.prior.create_default_prior(name='luminosity_distance')
priors['geocent_time'] = tupak.core.prior.Uniform(injection_parameters['geocent_time'] - 1,
injection_parameters['geocent_time'] + 1,
'geocent_time')
diff --git a/examples/injection_examples/basic_tutorial_dist_time_phase_marg.py b/examples/injection_examples/basic_tutorial_dist_time_phase_marg.py
index bad5ae9a5..0577b3828 100644
--- a/examples/injection_examples/basic_tutorial_dist_time_phase_marg.py
+++ b/examples/injection_examples/basic_tutorial_dist_time_phase_marg.py
@@ -54,7 +54,7 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl','psi', 'ra', 'd
# The above list does *not* include mass_1, mass_2, iota and luminosity_distance, which means those are the parameters
# that will be included in the sampler. If we do nothing, then the default priors get used.
-priors['luminosity_distance'] = tupak.gw.prior.create_default_prior(name='luminosity_distance')
+priors['luminosity_distance'] = tupak.core.prior.create_default_prior(name='luminosity_distance')
# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
likelihood = tupak.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,
diff --git a/examples/injection_examples/basic_tutorial_time_phase_marg.py b/examples/injection_examples/basic_tutorial_time_phase_marg.py
index 81c2670a9..a601d6702 100644
--- a/examples/injection_examples/basic_tutorial_time_phase_marg.py
+++ b/examples/injection_examples/basic_tutorial_time_phase_marg.py
@@ -54,7 +54,7 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl','psi', 'ra', 'd
# The above list does *not* include mass_1, mass_2, iota and luminosity_distance, which means those are the parameters
# that will be included in the sampler. If we do nothing, then the default priors get used.
-priors['luminosity_distance'] = tupak.gw.prior.create_default_prior(name='luminosity_distance')
+priors['luminosity_distance'] = tupak.core.prior.create_default_prior(name='luminosity_distance')
# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
likelihood = tupak.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,
diff --git a/examples/injection_examples/change_sampled_parameters.py b/examples/injection_examples/change_sampled_parameters.py
index 65f818d30..5ef99386e 100644
--- a/examples/injection_examples/change_sampled_parameters.py
+++ b/examples/injection_examples/change_sampled_parameters.py
@@ -41,7 +41,7 @@ priors = dict()
# These parameters will not be sampled
for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'psi', 'ra', 'dec', 'geocent_time']:
priors[key] = injection_parameters[key]
-priors['luminosity_distance'] = tupak.gw.prior.create_default_prior(name='luminosity_distance')
+priors['luminosity_distance'] = tupak.core.prior.create_default_prior(name='luminosity_distance')
# Initialise GravitationalWaveTransient
likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator)
diff --git a/examples/injection_examples/how_to_specify_the_prior.py b/examples/injection_examples/how_to_specify_the_prior.py
index 2c14b1e5d..e6d3aa8c1 100644
--- a/examples/injection_examples/how_to_specify_the_prior.py
+++ b/examples/injection_examples/how_to_specify_the_prior.py
@@ -38,7 +38,7 @@ priors = dict()
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 assign a default prior distribution, note this only works for certain parameters.
-priors['mass_1'] = tupak.gw.prior.create_default_prior(name='mass_1')
+priors['mass_1'] = tupak.core.prior.create_default_prior(name='mass_1')
# We can make uniform distributions.
priors['mass_2'] = tupak.core.prior.Uniform(name='mass_2', minimum=20, maximum=40)
# We can load a prior distribution from a file, e.g., a uniform in comoving volume distribution.
diff --git a/examples/supernova_example/supernova_example.py b/examples/supernova_example/supernova_example.py
index 1a194a965..9bfcbfe39 100644
--- a/examples/supernova_example/supernova_example.py
+++ b/examples/supernova_example/supernova_example.py
@@ -75,8 +75,8 @@ priors['pc_coeff2'] = tupak.core.prior.Uniform(-1, 1, 'pc_coeff2')
priors['pc_coeff3'] = tupak.core.prior.Uniform(-1, 1, 'pc_coeff3')
priors['pc_coeff4'] = tupak.core.prior.Uniform(-1, 1, 'pc_coeff4')
priors['pc_coeff5'] = tupak.core.prior.Uniform(-1, 1, 'pc_coeff5')
-priors['ra'] = tupak.gw.prior.create_default_prior(name='ra')
-priors['dec'] = tupak.gw.prior.create_default_prior(name='dec')
+priors['ra'] = tupak.core.prior.create_default_prior(name='ra')
+priors['dec'] = tupak.core.prior.create_default_prior(name='dec')
priors['geocent_time'] = tupak.core.prior.Uniform(
injection_parameters['geocent_time'] - 1,
injection_parameters['geocent_time'] + 1,
diff --git a/tupak/core/prior.py b/tupak/core/prior.py
index 4a9fb470e..811324af9 100644
--- a/tupak/core/prior.py
+++ b/tupak/core/prior.py
@@ -9,7 +9,7 @@ import scipy.stats
import logging
import os
-from tupak.gw.prior import create_default_prior, test_redundancy
+from tupak.gw.prior import test_redundancy
class PriorSet(dict):
@@ -62,6 +62,96 @@ class PriorSet(dict):
prior[key] = eval(val)
self.update(prior)
+ def fill_priors(self, likelihood):
+ """
+ Fill dictionary of priors based on required parameters of likelihood
+
+ Any floats in prior will be converted to delta function prior. Any
+ required, non-specified parameters will use the default.
+
+ Parameters
+ ----------
+ prior: dict
+ dictionary of prior objects and floats
+ likelihood: tupak.likelihood.GravitationalWaveTransient instance
+ Used to infer the set of parameters to fill the prior with
+
+ Note: if `likelihood` has `non_standard_sampling_parameter_keys`, then
+ this will set-up default priors for those as well.
+
+ Returns
+ -------
+ prior: dict
+ The filled prior dictionary
+
+ """
+
+ for key in self:
+ if isinstance(self[key], Prior):
+ continue
+ elif isinstance(self[key], float) or isinstance(self[key], int):
+ self[key] = DeltaFunction(self[key])
+ logging.info(
+ "{} converted to delta function prior.".format(key))
+ else:
+ logging.info(
+ "{} cannot be converted to delta function prior."
+ .format(key))
+
+ missing_keys = set(likelihood.parameters) - set(self.keys())
+
+ if getattr(likelihood, 'non_standard_sampling_parameter_keys', None) is not None:
+ for parameter in likelihood.non_standard_sampling_parameter_keys:
+ self[parameter] = create_default_prior(parameter)
+
+ for missing_key in missing_keys:
+ default_prior = create_default_prior(missing_key)
+ if default_prior is None:
+ set_val = likelihood.parameters[missing_key]
+ logging.warning(
+ "Parameter {} has no default prior and is set to {}, this"
+ " will not be sampled and may cause an error."
+ .format(missing_key, set_val))
+ else:
+ if not test_redundancy(missing_key, self):
+ self[missing_key] = default_prior
+
+ for key in self:
+ test_redundancy(key, self)
+
+
+def create_default_prior(name, default_prior_file=None):
+ """
+ Make a default prior for a parameter with a known name.
+
+ Parameters
+ ----------
+ name: str
+ Parameter name
+ default_prior_file: str
+ If given, the file where default priors are stored.
+
+ Return
+ ------
+ prior: Prior
+ Default prior distribution for that parameter, if unknown None is
+ returned.
+ """
+
+ if default_prior_file is None:
+ default_prior_file = 'binary_black_holes.prior'
+ default_prior_file = os.path.join(os.path.dirname(__file__),
+ 'prior_files',
+ 'binary_black_holes.prior')
+ default_priors = PriorSet(filename=default_prior_file)
+ if name in default_priors.keys():
+ prior = default_priors[name]
+ else:
+ logging.info(
+ "No default prior found for variable {}.".format(name))
+ prior = None
+ return prior
+
class Prior(object):
"""
@@ -501,61 +591,4 @@ class UniformComovingVolume(FromFile):
return FromFile.__repr__(self)
-def fill_priors(prior, likelihood):
- """
- Fill dictionary of priors based on required parameters of likelihood
-
- Any floats in prior will be converted to delta function prior. Any
- required, non-specified parameters will use the default.
-
- Parameters
- ----------
- prior: dict
- dictionary of prior objects and floats
- likelihood: tupak.likelihood.GravitationalWaveTransient instance
- Used to infer the set of parameters to fill the prior with
-
- Note: if `likelihood` has `non_standard_sampling_parameter_keys`, then this
- will set-up default priors for those as well.
-
- Returns
- -------
- prior: dict
- The filled prior dictionary
-
- """
-
- for key in prior:
- if isinstance(prior[key], Prior):
- continue
- elif isinstance(prior[key], float) or isinstance(prior[key], int):
- prior[key] = DeltaFunction(prior[key])
- logging.info(
- "{} converted to delta function prior.".format(key))
- else:
- logging.info(
- "{} cannot be converted to delta function prior.".format(key))
-
- missing_keys = set(likelihood.parameters) - set(prior.keys())
-
- if getattr(likelihood, 'non_standard_sampling_parameter_keys', None) is not None:
- for parameter in likelihood.non_standard_sampling_parameter_keys:
- prior[parameter] = create_default_prior(parameter)
-
- for missing_key in missing_keys:
- default_prior = create_default_prior(missing_key)
- if default_prior is None:
- set_val = likelihood.parameters[missing_key]
- logging.warning(
- "Parameter {} has no default prior and is set to {}, this will"
- " not be sampled and may cause an error."
- .format(missing_key, set_val))
- else:
- if not test_redundancy(missing_key, prior):
- prior[missing_key] = default_prior
-
- for key in prior:
- test_redundancy(key, prior)
-
- return prior
diff --git a/tupak/core/prior_files/binary_black_holes.prior b/tupak/core/prior_files/binary_black_holes.prior
new file mode 100644
index 000000000..ccb699d64
--- /dev/null
+++ b/tupak/core/prior_files/binary_black_holes.prior
@@ -0,0 +1,21 @@
+mass_1 = Uniform(name='mass_1', minimum=20, maximum=100)
+mass_2 = Uniform(name='mass_2', minimum=20, maximum=100)
+chirp_mass = Uniform(name='chirp_mass', minimum=25, maximum=100)
+total_mass = Uniform(name='total_mass', minimum=10, maximum=200)
+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)
+a_2 = Uniform(name='a_2', minimum=0, maximum=0.8)
+tilt_1 = Sine(name='tilt_1')
+tilt_2 = Sine(name='tilt_2')
+cos_tilt_1 = Uniform(name='cos_tilt_1', minimum=-1, maximum=1)
+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 = UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3)
+dec = Cosine(name='dec')
+ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi)
+iota = Sine(name='iota')
+cos_iota = Uniform(name='cos_iota', minimum=-1, maximum=1)
+psi = Uniform(name='psi', minimum=0, maximum=2 * np.pi)
+phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi)
diff --git a/tupak/core/sampler.py b/tupak/core/sampler.py
index 6c4d24676..1fcaf5ee5 100644
--- a/tupak/core/sampler.py
+++ b/tupak/core/sampler.py
@@ -9,7 +9,7 @@ import matplotlib.pyplot as plt
import time
from tupak.core.result import Result, read_in_result
-from tupak.core.prior import Prior, fill_priors
+from tupak.core.prior import Prior
from tupak.core import utils
import tupak
@@ -556,7 +556,6 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
if priors is None:
priors = dict()
- priors = fill_priors(priors, likelihood)
if type(priors) == dict:
priors = tupak.core.prior.PriorSet(priors)
@@ -565,6 +564,7 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
else:
raise ValueError
+ priors.fill_priors(likelihood)
priors.write_to_file(outdir, label)
if implemented_samplers.__contains__(sampler.title()):
diff --git a/tupak/gw/likelihood.py b/tupak/gw/likelihood.py
index d01334468..cee3ce703 100644
--- a/tupak/gw/likelihood.py
+++ b/tupak/gw/likelihood.py
@@ -178,7 +178,7 @@ class GravitationalWaveTransient(likelihood.Likelihood):
def setup_distance_marginalization(self):
if 'luminosity_distance' not in self.prior.keys():
logging.info('No prior provided for distance, using default prior.')
- self.prior['luminosity_distance'] = tupak.gw.prior.create_default_prior('luminosity_distance')
+ self.prior['luminosity_distance'] = tupak.core.prior.create_default_prior('luminosity_distance')
self.distance_array = np.linspace(self.prior['luminosity_distance'].minimum,
self.prior['luminosity_distance'].maximum, int(1e4))
self.delta_distance = self.distance_array[1] - self.distance_array[0]
@@ -219,7 +219,7 @@ class GravitationalWaveTransient(likelihood.Likelihood):
def setup_phase_marginalization(self):
if 'phase' not in self.prior.keys() or not isinstance(self.prior['phase'], tupak.core.prior.Prior):
logging.info('No prior provided for phase at coalescence, using default prior.')
- self.prior['phase'] = tupak.gw.prior.create_default_prior('phase')
+ self.prior['phase'] = tupak.core.prior.create_default_prior('phase')
self.bessel_function_interped = interp1d(np.linspace(0, 1e6, int(1e5)),
np.log([i0e(snr) for snr in np.linspace(0, 1e6, int(1e5))])
+ np.linspace(0, 1e6, int(1e5)),
diff --git a/tupak/gw/prior.py b/tupak/gw/prior.py
index cabe97bea..456a6ab7a 100644
--- a/tupak/gw/prior.py
+++ b/tupak/gw/prior.py
@@ -1,58 +1,8 @@
import logging
-import numpy as np
-
import tupak.core.prior
-def create_default_prior(name):
- """
- Make a default prior for a parameter with a known name.
-
- This is currently set up for binary black holes.
-
- Parameters
- ----------
- name: str
- Parameter name
-
- Return
- ------
- prior: Prior
- Default prior distribution for that parameter, if unknown None is returned.
- """
- default_priors = {
- 'mass_1': tupak.core.prior.Uniform(name=name, minimum=20, maximum=100),
- 'mass_2': tupak.core.prior.Uniform(name=name, minimum=20, maximum=100),
- 'chirp_mass': tupak.core.prior.Uniform(name=name, minimum=25, maximum=100),
- 'total_mass': tupak.core.prior.Uniform(name=name, minimum=10, maximum=200),
- 'mass_ratio': tupak.core.prior.Uniform(name=name, minimum=0.125, maximum=1),
- 'symmetric_mass_ratio': tupak.core.prior.Uniform(name=name, minimum=8 / 81, maximum=0.25),
- 'a_1': tupak.core.prior.Uniform(name=name, minimum=0, maximum=0.8),
- 'a_2': tupak.core.prior.Uniform(name=name, minimum=0, maximum=0.8),
- 'tilt_1': tupak.core.prior.Sine(name=name),
- 'tilt_2': tupak.core.prior.Sine(name=name),
- 'cos_tilt_1': tupak.core.prior.Uniform(name=name, minimum=-1, maximum=1),
- 'cos_tilt_2': tupak.core.prior.Uniform(name=name, minimum=-1, maximum=1),
- 'phi_12': tupak.core.prior.Uniform(name=name, minimum=0, maximum=2 * np.pi),
- 'phi_jl': tupak.core.prior.Uniform(name=name, minimum=0, maximum=2 * np.pi),
- 'luminosity_distance': tupak.core.prior.UniformComovingVolume(name=name, minimum=1e2, maximum=5e3),
- 'dec': tupak.core.prior.Cosine(name=name),
- 'ra': tupak.core.prior.Uniform(name=name, minimum=0, maximum=2 * np.pi),
- 'iota': tupak.core.prior.Sine(name=name),
- 'cos_iota': tupak.core.prior.Uniform(name=name, minimum=-1, maximum=1),
- 'psi': tupak.core.prior.Uniform(name=name, minimum=0, maximum=2 * np.pi),
- 'phase': tupak.core.prior.Uniform(name=name, minimum=0, maximum=2 * np.pi)
- }
- if name in default_priors.keys():
- prior = default_priors[name]
- else:
- logging.info(
- "No default prior found for variable {}.".format(name))
- prior = None
- return prior
-
-
def test_redundancy(key, prior):
"""
Test whether adding the key would add be redundant.
@@ -99,4 +49,4 @@ def test_redundancy(key, prior):
redundant = True
break
- return redundant
\ No newline at end of file
+ return redundant
--
GitLab