diff --git a/examples/injection_examples/basic_tutorial.py b/examples/injection_examples/basic_tutorial.py index 85cb46977fe00072ed3e95d5998f153d0039fd9f..00e8403b272594f1bde053058d1088a3394f107d 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 bad5ae9a542ea4e45d970fb4fcc5e817f91d261c..0577b382843790b1407d0420fc7031bcff901c3a 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 81c2670a927813d1f41312ba50fa168dd90572cb..a601d670261ad3917e2d9204b5ce2bcbf19d03ab 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 65f818d30be3087b7b905b9520f3bd9b95aca4b4..5ef99386ea412aebb9fa910d3a7ae58c5c3f247e 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 2c14b1e5d84c16560d623d7bce9d5083009b2236..e6d3aa8c12209ea2d4387e22e6e2c44b6831ba27 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 1a194a965bc3426eccf95ae6d253b332f16dde37..9bfcbfe39b5f46884c1f014bb533d96ed613f938 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 4a9fb470e9068574bf4f608698ef908039add9f9..811324af9a1ce26d0fc867b82c63f7960c2844f2 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 0000000000000000000000000000000000000000..ccb699d64a499d31cb9b7a8a7ceacdfb0bef524b --- /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 6c4d24676ae93a73b3cfeff75f2aea0ad5e08865..1fcaf5ee59403a42d6920214fbc3dedc92167e85 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 d01334468456639fb5a11c2d36f6430d180501c9..cee3ce703cd7bc7a75da764b76eb76d831680f58 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 cabe97beaf8eca0b7facb455a9b24050dfc835b4..456a6ab7a8ea3c137f877b67a7979a88928388a0 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