From cc13f3e8e564ebb32b9b1b85f65d744d7cf63ed8 Mon Sep 17 00:00:00 2001 From: Colm Talbot <colm.talbot@ligo.org> Date: Thu, 4 Oct 2018 21:14:42 +1000 Subject: [PATCH] pep examples --- examples/injection_examples/basic_tutorial.py | 64 +++++++++++-------- .../change_sampled_parameters.py | 6 ++ .../create_your_own_source_model.py | 16 +++-- ...reate_your_own_time_domain_source_model.py | 42 ++++++------ .../injection_examples/eccentric_inspiral.py | 45 ++++++------- .../how_to_specify_the_prior.py | 30 +++++---- .../marginalized_likelihood.py | 31 +++++---- .../sine_gaussian_example.py | 48 +++++++------- examples/injection_examples/using_gwin.py | 4 +- .../other_examples/linear_regression_pymc3.py | 4 +- ...near_regression_pymc3_custom_likelihood.py | 8 ++- .../supernova_example/supernova_example.py | 8 +-- 12 files changed, 167 insertions(+), 139 deletions(-) diff --git a/examples/injection_examples/basic_tutorial.py b/examples/injection_examples/basic_tutorial.py index 1e90bcddc..e22d38056 100644 --- a/examples/injection_examples/basic_tutorial.py +++ b/examples/injection_examples/basic_tutorial.py @@ -1,18 +1,19 @@ #!/usr/bin/env python """ -Tutorial to demonstrate running parameter estimation on a reduced parameter space for an injected signal. +Tutorial to demonstrate running parameter estimation on a reduced parameter +space for an injected signal. -This example estimates the masses using a uniform prior in both component masses and distance using a uniform in -comoving volume prior on luminosity distance between luminosity distances of 100Mpc and 5Gpc, the cosmology is WMAP7. +This example estimates the masses using a uniform prior in both component masses +and distance using a uniform in comoving volume prior on luminosity distance +between luminosity distances of 100Mpc and 5Gpc, the cosmology is WMAP7. """ from __future__ import division, print_function import numpy as np - import bilby -# Set the duration and sampling frequency of the data segment that we're going to inject the signal into - +# Set the duration and sampling frequency of the data segment that we're +# going to inject the signal into duration = 4. sampling_frequency = 2048. @@ -24,22 +25,24 @@ bilby.core.utils.setup_logger(outdir=outdir, label=label) # Set up a random seed for result reproducibility. This is optional! np.random.seed(88170235) -# We are going to inject a binary black hole waveform. We first establish a dictionary of parameters that -# includes all of the different waveform parameters, including masses of the two black holes (mass_1, mass_2), +# We are going to inject a binary black hole waveform. We first establish a +# dictionary of parameters that includes all of the different waveform +# parameters, including masses of the two black holes (mass_1, mass_2), # spins of both black holes (a, tilt, phi), etc. -injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, phi_12=1.7, phi_jl=0.3, - luminosity_distance=2000., iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413, - ra=1.375, dec=-1.2108) +injection_parameters = dict( + mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, + phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., iota=0.4, psi=2.659, + phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108) # Fixed arguments passed into the source model waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', - reference_frequency=50.) + reference_frequency=50., minimum_frequency=20.) # Create the waveform_generator using a LAL BinaryBlackHole source function waveform_generator = bilby.gw.WaveformGenerator( duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole, - parameters=injection_parameters, waveform_arguments=waveform_arguments) + waveform_arguments=waveform_arguments) # Set up interferometers. In this case we'll use two interferometers # (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design @@ -51,30 +54,35 @@ ifos.set_strain_data_from_power_spectral_densities( ifos.inject_signal(waveform_generator=waveform_generator, parameters=injection_parameters) -# Set up prior, which is a dictionary -# By default we will sample all terms in the signal models. However, this will take a long time for the calculation, -# so for this example we will set almost all of the priors to be equall to their injected values. This implies the -# prior is a delta function at the true, injected value. In reality, the sampler implementation is smart enough to -# not sample any parameter that has a delta-function prior. -# 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. +# Set up a PriorSet, which inherits from dict. +# By default we will sample all terms in the signal models. However, this will +# take a long time for the calculation, so for this example we will set almost +# all of the priors to be equall to their injected values. This implies the +# prior is a delta function at the true, injected value. In reality, the +# sampler implementation is smart enough to not sample any parameter that has +# a delta-function prior. +# 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 = bilby.gw.prior.BBHPriorSet() priors['geocent_time'] = bilby.core.prior.Uniform( 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']: +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] -# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator -likelihood = bilby.gw.GravitationalWaveTransient(interferometers=ifos, waveform_generator=waveform_generator, - time_marginalization=False, phase_marginalization=False, - distance_marginalization=False, prior=priors) +# Initialise the likelihood by passing in the interferometer data (ifos) and +# the waveoform generator +likelihood = bilby.gw.GravitationalWaveTransient( + interferometers=ifos, waveform_generator=waveform_generator) # Run sampler. In this case we're going to use the `dynesty` sampler -result = bilby.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000, - injection_parameters=injection_parameters, outdir=outdir, label=label) +result = bilby.run_sampler( + likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000, + injection_parameters=injection_parameters, outdir=outdir, label=label) -# make some plots of the outputs +# Make a corner plot. result.plot_corner() diff --git a/examples/injection_examples/change_sampled_parameters.py b/examples/injection_examples/change_sampled_parameters.py index f629b2c91..aabcbfebc 100644 --- a/examples/injection_examples/change_sampled_parameters.py +++ b/examples/injection_examples/change_sampled_parameters.py @@ -29,6 +29,8 @@ waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', reference_frequency=50.) # Create the waveform_generator using a LAL BinaryBlackHole source function +# We specify a function which transforms a dictionary of parameters into the +# appropriate parameters for the source model. waveform_generator = bilby.gw.waveform_generator.WaveformGenerator( sampling_frequency=sampling_frequency, duration=duration, frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole, @@ -45,6 +47,8 @@ ifos.inject_signal(waveform_generator=waveform_generator, parameters=injection_parameters) # Set up prior +# Note it is possible to sample in different parameters to those that were +# injected. priors = bilby.gw.prior.BBHPriorSet() priors.pop('mass_1') priors.pop('mass_2') @@ -69,6 +73,8 @@ likelihood = bilby.gw.likelihood.GravitationalWaveTransient( interferometers=ifos, waveform_generator=waveform_generator) # Run sampler +# Note we've added a post-processing conversion function, this will generate +# many useful additional parameters, e.g., source-frame masses. result = bilby.core.sampler.run_sampler( likelihood=likelihood, priors=priors, sampler='dynesty', outdir=outdir, injection_parameters=injection_parameters, label='DifferentParameters', diff --git a/examples/injection_examples/create_your_own_source_model.py b/examples/injection_examples/create_your_own_source_model.py index 254fdfc52..4e0ab71da 100644 --- a/examples/injection_examples/create_your_own_source_model.py +++ b/examples/injection_examples/create_your_own_source_model.py @@ -22,13 +22,14 @@ def sine_gaussian(f, A, f0, tau, phi0, geocent_time, ra, dec, psi): return {'plus': plus, 'cross': cross} -# We now define some parameters that we will inject and then a waveform generator +# We now define some parameters that we will inject injection_parameters = dict(A=1e-23, f0=100, tau=1, phi0=0, geocent_time=0, ra=0, dec=0, psi=0) -waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(duration=duration, - sampling_frequency=sampling_frequency, - frequency_domain_source_model=sine_gaussian, - parameters=injection_parameters) + +# Now we pass our source function to the WaveformGenerator +waveform_generator = bilby.gw.waveform_generator.WaveformGenerator( + duration=duration, sampling_frequency=sampling_frequency, + frequency_domain_source_model=sine_gaussian) # Set up interferometers. ifos = bilby.gw.detector.InterferometerList(['H1', 'L1']) @@ -41,10 +42,11 @@ ifos.inject_signal(waveform_generator=waveform_generator, # Here we define the priors for the search. We use the injection parameters # except for the amplitude, f0, and geocent_time prior = injection_parameters.copy() -prior['A'] = bilby.core.prior.PowerLaw(alpha=-1, minimum=1e-25, maximum=1e-21, name='A') +prior['A'] = bilby.core.prior.LogUniform(minimum=1e-25, maximum=1e-21, name='A') prior['f0'] = bilby.core.prior.Uniform(90, 110, 'f') -likelihood = bilby.gw.likelihood.GravitationalWaveTransient(ifos, waveform_generator) +likelihood = bilby.gw.likelihood.GravitationalWaveTransient( + interferometers=ifos, waveform_generator=waveform_generator) result = bilby.core.sampler.run_sampler( likelihood, prior, sampler='dynesty', outdir=outdir, label=label, diff --git a/examples/injection_examples/create_your_own_time_domain_source_model.py b/examples/injection_examples/create_your_own_time_domain_source_model.py index 5ea57fb58..28bae6952 100644 --- a/examples/injection_examples/create_your_own_time_domain_source_model.py +++ b/examples/injection_examples/create_your_own_time_domain_source_model.py @@ -1,41 +1,41 @@ #!/usr/bin/env python """ A script to show how to create your own time domain source model. -A simple damped Gaussian signal is defined in the time domain, injected into noise in -two interferometers (LIGO Livingston and Hanford at design sensitivity), -and then recovered. +A simple damped Gaussian signal is defined in the time domain, injected into +noise in two interferometers (LIGO Livingston and Hanford at design +sensitivity), and then recovered. """ -import bilby import numpy as np - +import bilby # define the time-domain model -def time_domain_damped_sinusoid(time, amplitude, damping_time, frequency, phase, ra, dec, psi, geocent_time): +def time_domain_damped_sinusoid( + time, amplitude, damping_time, frequency, phase): """ - This example only creates a linearly polarised signal with only plus polarisation. + This example only creates a linearly polarised signal with only plus + polarisation. """ - - plus = amplitude * np.exp(-time / damping_time) * np.sin(2.*np.pi*frequency*time + phase) + plus = amplitude * np.exp(-time / damping_time) *\ + np.sin(2 * np.pi * frequency * time + phase) cross = np.zeros(len(time)) - return {'plus': plus, 'cross': cross} + # define parameters to inject. injection_parameters = dict(amplitude=5e-22, damping_time=0.1, frequency=50, - phase=0, - ra=0, dec=0, psi=0, geocent_time=0.) + phase=0, ra=0, dec=0, psi=0, geocent_time=0.) duration = 0.5 sampling_frequency = 2048 -outdir='outdir' -label='time_domain_source_model' +outdir = 'outdir' +label = 'time_domain_source_model' # call the waveform_generator to create our waveform model. -waveform = bilby.gw.waveform_generator.WaveformGenerator(duration=duration, sampling_frequency=sampling_frequency, - time_domain_source_model=time_domain_damped_sinusoid, - parameters=injection_parameters) +waveform = bilby.gw.waveform_generator.WaveformGenerator( + duration=duration, sampling_frequency=sampling_frequency, + time_domain_source_model=time_domain_damped_sinusoid) # inject the signal into three interferometers ifos = bilby.gw.detector.InterferometerList(['H1', 'L1']) @@ -47,7 +47,7 @@ ifos.inject_signal(waveform_generator=waveform, # create the priors prior = injection_parameters.copy() -prior['amplitude'] = bilby.core.prior.Uniform(1e-23, 1e-21, r'$h_0$') +prior['amplitude'] = bilby.core.prior.LogUniform(1e-23, 1e-21, r'$h_0$') prior['damping_time'] = bilby.core.prior.Uniform( 0, 1, r'damping time', unit='$s$') prior['frequency'] = bilby.core.prior.Uniform(0, 200, r'frequency', unit='Hz') @@ -57,9 +57,9 @@ prior['phase'] = bilby.core.prior.Uniform(-np.pi / 2, np.pi / 2, r'$\phi$') likelihood = bilby.gw.likelihood.GravitationalWaveTransient(ifos, waveform) # launch sampler -result = bilby.core.sampler.run_sampler(likelihood, prior, sampler='dynesty', npoints=1000, - injection_parameters=injection_parameters, - outdir=outdir, label=label) +result = bilby.core.sampler.run_sampler( + likelihood, prior, sampler='dynesty', npoints=1000, + injection_parameters=injection_parameters, outdir=outdir, label=label) result.plot_corner() diff --git a/examples/injection_examples/eccentric_inspiral.py b/examples/injection_examples/eccentric_inspiral.py index 75ff93b57..2034df2b5 100644 --- a/examples/injection_examples/eccentric_inspiral.py +++ b/examples/injection_examples/eccentric_inspiral.py @@ -1,8 +1,8 @@ #!/usr/bin/env python """ -Tutorial to demonstrate running parameter estimation on a reduced parameter space -for an injected eccentric binary black hole signal with masses & distnace similar -to GW150914. +Tutorial to demonstrate running parameter estimation on a reduced parameter +space for an injected eccentric binary black hole signal with masses & distnace +similar to GW150914. This uses the same binary parameters that were used to make Figures 1, 2 & 5 in Lower et al. (2018) -> arXiv:1806.05350. @@ -10,14 +10,11 @@ Lower et al. (2018) -> arXiv:1806.05350. For a more comprehensive look at what goes on in each step, refer to the "basic_tutorial.py" example. """ -from __future__ import division, print_function +from __future__ import division import numpy as np - import bilby -import matplotlib.pyplot as plt - duration = 64. sampling_frequency = 256. @@ -28,13 +25,15 @@ bilby.core.utils.setup_logger(outdir=outdir, label=label) # Set up a random seed for result reproducibility. np.random.seed(150914) -injection_parameters = dict(mass_1=35., mass_2=30., eccentricity=0.1, - luminosity_distance=440., iota=0.4, psi=0.1, phase=1.2, - geocent_time=1180002601.0, ra=45, dec=5.73) +injection_parameters = dict( + mass_1=35., mass_2=30., eccentricity=0.1, luminosity_distance=440., + iota=0.4, psi=0.1, phase=1.2, geocent_time=1180002601.0, ra=45, dec=5.73) -waveform_arguments = dict(waveform_approximant='EccentricFD', reference_frequency=10., minimum_frequency=10.) +waveform_arguments = dict(waveform_approximant='EccentricFD', + reference_frequency=10., minimum_frequency=10.) -# Create the waveform_generator using the LAL eccentric black hole no spins source function +# Create the waveform_generator using the LAL eccentric black hole no spins +# source function waveform_generator = bilby.gw.WaveformGenerator( duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=bilby.gw.source.lal_eccentric_binary_black_hole_no_spins, @@ -43,8 +42,8 @@ waveform_generator = bilby.gw.WaveformGenerator( # Setting up three interferometers (LIGO-Hanford (H1), LIGO-Livingston (L1), and # Virgo (V1)) at their design sensitivities. The maximum frequency is set just -# prior to the point at which the waveform model terminates. This is to avoid any -# biases introduced from using a sharply terminating waveform model. +# prior to the point at which the waveform model terminates. This is to avoid +# any biases introduced from using a sharply terminating waveform model. minimum_frequency = 10. maximum_frequency = 128. @@ -59,13 +58,13 @@ ifos.inject_signal(waveform_generator=waveform_generator, parameters=injection_parameters) # Now we set up the priors on each of the binary parameters. -priors = dict() +priors = bilby.core.prior.PriorSet() priors["mass_1"] = bilby.core.prior.Uniform( name='mass_1', minimum=5, maximum=60, unit='$M_{\\odot}$') priors["mass_2"] = bilby.core.prior.Uniform( name='mass_2', minimum=5, maximum=60, unit='$M_{\\odot}$') -priors["eccentricity"] = bilby.core.prior.PowerLaw( - name='eccentricity', latex_label='$e$', alpha=-1, minimum=1e-4, maximum=0.4) +priors["eccentricity"] = bilby.core.prior.LogUniform( + name='eccentricity', latex_label='$e$', minimum=1e-4, maximum=0.4) priors["luminosity_distance"] = bilby.gw.prior.UniformComovingVolume( name='luminosity_distance', minimum=1e2, maximum=2e3) priors["dec"] = bilby.core.prior.Cosine(name='dec') @@ -79,15 +78,13 @@ priors["geocent_time"] = bilby.core.prior.Uniform( 1180002600.9, 1180002601.1, name='geocent_time', unit='s') # Initialising the likelihood function. -likelihood = bilby.gw.likelihood.GravitationalWaveTransient(interferometers=ifos, - waveform_generator=waveform_generator, time_marginalization=False, - phase_marginalization=False, distance_marginalization=False, - prior=priors) +likelihood = bilby.gw.likelihood.GravitationalWaveTransient( + interferometers=ifos, waveform_generator=waveform_generator) # Now we run sampler (PyMultiNest in our case). -result = bilby.run_sampler(likelihood=likelihood, priors=priors, sampler='pymultinest', - npoints=1000, injection_parameters=injection_parameters, - outdir=outdir, label=label) +result = bilby.run_sampler( + likelihood=likelihood, priors=priors, sampler='pymultinest', npoints=1000, + injection_parameters=injection_parameters, outdir=outdir, label=label) # And finally we make some plots of the output posteriors. result.plot_corner() diff --git a/examples/injection_examples/how_to_specify_the_prior.py b/examples/injection_examples/how_to_specify_the_prior.py index b92288d90..4c9b80635 100644 --- a/examples/injection_examples/how_to_specify_the_prior.py +++ b/examples/injection_examples/how_to_specify_the_prior.py @@ -1,12 +1,12 @@ #!/usr/bin/env python """ -Tutorial to demonstrate how to specify the prior distributions used for parameter estimation. +Tutorial to demonstrate how to specify the prior distributions used for +parameter estimation. """ from __future__ import division, print_function -import bilby -import numpy as np -import bilby.gw.prior +import numpy as np +import bilby duration = 4. @@ -15,18 +15,19 @@ outdir = 'outdir' np.random.seed(151012) -injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, phi_12=1.7, phi_jl=0.3, - luminosity_distance=4000., iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413, - ra=1.375, dec=-1.2108) +injection_parameters = dict( + mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, + phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4, psi=2.659, + phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108) waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', - reference_frequency=50.) + reference_frequency=50., minimum_frequency=20.) # Create the waveform_generator using a LAL BinaryBlackHole source function waveform_generator = bilby.gw.WaveformGenerator( duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole, - parameters=injection_parameters, waveform_arguments=waveform_arguments) + waveform_arguments=waveform_arguments) # Set up interferometers. ifos = bilby.gw.detector.InterferometerList(['H1', 'L1']) @@ -58,14 +59,17 @@ priors['a_2'] = bilby.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. +# 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 = bilby.gw.GravitationalWaveTransient(interferometers=ifos, waveform_generator=waveform_generator) +likelihood = bilby.gw.GravitationalWaveTransient( + interferometers=ifos, waveform_generator=waveform_generator) # Run sampler -result = bilby.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', - injection_parameters=injection_parameters, outdir=outdir, label='specify_prior') +result = bilby.run_sampler( + likelihood=likelihood, priors=priors, sampler='dynesty', outdir=outdir, + injection_parameters=injection_parameters, label='specify_prior') result.plot_corner() diff --git a/examples/injection_examples/marginalized_likelihood.py b/examples/injection_examples/marginalized_likelihood.py index 83259d627..29eae13c6 100644 --- a/examples/injection_examples/marginalized_likelihood.py +++ b/examples/injection_examples/marginalized_likelihood.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ -Tutorial to demonstrate how to improve the speed and efficiency of parameter estimation on an injected signal using -phase and distance marginalisation. +Tutorial to demonstrate how to improve the speed and efficiency of parameter +estimation on an injected signal using time, phase and distance marginalisation. """ from __future__ import division, print_function import bilby @@ -11,12 +11,14 @@ import numpy as np duration = 4. sampling_frequency = 2048. outdir = 'outdir' +label = 'marginalized_likelihood' np.random.seed(170608) -injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, phi_12=1.7, phi_jl=0.3, - luminosity_distance=4000., iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413, - ra=1.375, dec=-1.2108) +injection_parameters = dict( + mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, + phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4, psi=2.659, + phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108) waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', reference_frequency=50.) @@ -24,7 +26,7 @@ waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', # Create the waveform_generator using a LAL BinaryBlackHole source function waveform_generator = bilby.gw.WaveformGenerator( duration=duration, sampling_frequency=sampling_frequency, - frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole, parameters=injection_parameters, + frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole, waveform_arguments=waveform_arguments) # Set up interferometers. @@ -38,19 +40,22 @@ ifos.inject_signal(waveform_generator=waveform_generator, # Set up prior priors = bilby.gw.prior.BBHPriorSet() # These parameters will not be sampled -for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'iota', 'ra', 'dec', 'geocent_time']: +for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'iota', 'ra', + 'dec']: priors[key] = injection_parameters[key] # Initialise GravitationalWaveTransient -# Note that we now need to pass the: priors and flags for each thing that's being marginalised. -# This is still under development so care should be taken with the marginalised likelihood. +# Note that we now need to pass the: priors and flags for each thing that's +# being marginalised. A lookup table is used fro distance marginalisation which +# takes a few minutes to build. likelihood = bilby.gw.GravitationalWaveTransient( interferometers=ifos, waveform_generator=waveform_generator, prior=priors, - distance_marginalization=False, phase_marginalization=True, - time_marginalization=False) + distance_marginalization=True, phase_marginalization=True, + time_marginalization=True) # Run sampler -result = bilby.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', - injection_parameters=injection_parameters, outdir=outdir, label='MarginalisedLikelihood') +result = bilby.run_sampler( + likelihood=likelihood, priors=priors, sampler='dynesty', + injection_parameters=injection_parameters, outdir=outdir, label=label) result.plot_corner() diff --git a/examples/injection_examples/sine_gaussian_example.py b/examples/injection_examples/sine_gaussian_example.py index 4bbf7c6b0..2f0d1a127 100644 --- a/examples/injection_examples/sine_gaussian_example.py +++ b/examples/injection_examples/sine_gaussian_example.py @@ -1,13 +1,14 @@ #!/usr/bin/env python """ -Tutorial to demonstrate running parameter estimation on a sine gaussian injected signal. - +Tutorial to demonstrate running parameter estimation on a sine gaussian +injected signal. """ from __future__ import division, print_function import bilby import numpy as np -# Set the duration and sampling frequency of the data segment that we're going to inject the signal into +# Set the duration and sampling frequency of the data segment that we're going +# to inject the signal into duration = 4. sampling_frequency = 2048. @@ -19,19 +20,21 @@ bilby.core.utils.setup_logger(outdir=outdir, label=label) # Set up a random seed for result reproducibility. This is optional! np.random.seed(170801) -# We are going to inject a sine gaussian waveform. We first establish a dictionary of parameters that -# includes all of the different waveform parameters -injection_parameters = dict(hrss = 1e-22, Q = 5.0, frequency = 200.0, ra = 1.375, dec = -1.2108, - geocent_time = 1126259642.413, psi= 2.659) +# We are going to inject a sine gaussian waveform. We first establish a +# dictionary of parameters that includes all of the different waveform +# parameters +injection_parameters = dict( + hrss=1e-22, Q=5.0, frequency=200.0, ra=1.375, dec=-1.2108, + geocent_time=1126259642.413, psi=2.659) # Create the waveform_generator using a sine Gaussian source function -waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(duration=duration, - sampling_frequency=sampling_frequency, - frequency_domain_source_model=bilby.gw.source.sinegaussian, - parameters=injection_parameters) +waveform_generator = bilby.gw.waveform_generator.WaveformGenerator( + duration=duration, sampling_frequency=sampling_frequency, + frequency_domain_source_model=bilby.gw.source.sinegaussian) -# Set up interferometers. In this case we'll use three interferometers (LIGO-Hanford (H1), LIGO-Livingston (L1), -# and Virgo (V1)). These default to their design sensitivity +# Set up interferometers. In this case we'll use three interferometers +# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)). These default to +# their design sensitivity ifos = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1']) ifos.set_strain_data_from_power_spectral_densities( sampling_frequency=sampling_frequency, duration=duration, @@ -41,27 +44,22 @@ ifos.inject_signal(waveform_generator=waveform_generator, # Set up prior, which is a dictionary priors = dict() -# By default we will sample all terms in the signal models. However, this will take a long time for the calculation, -# so for this example we will set almost all of the priors to be equall to their injected values. This implies the -# prior is a delta function at the true, injected value. In reality, the sampler implementation is smart enough to -# not sample any parameter that has a delta-function prior. for key in ['psi', 'ra', 'dec', 'geocent_time']: priors[key] = injection_parameters[key] -# The above list does *not* include frequency and Q, which means those are the parameters -# that will be included in the sampler. If we do nothing, then the default priors get used. -#priors['Q'] = bilby.prior.create_default_prior(name='Q') -#priors['frequency'] = bilby.prior.create_default_prior(name='frequency') priors['Q'] = bilby.core.prior.Uniform(2, 50, 'Q') priors['frequency'] = bilby.core.prior.Uniform(30, 1000, 'frequency', unit='Hz') priors['hrss'] = bilby.core.prior.Uniform(1e-23, 1e-21, 'hrss') -# Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator -likelihood = bilby.gw.likelihood.GravitationalWaveTransient(interferometers=ifos, waveform_generator=waveform_generator) +# Initialise the likelihood by passing in the interferometer data (IFOs) and +# the waveoform generator +likelihood = bilby.gw.likelihood.GravitationalWaveTransient( + interferometers=ifos, waveform_generator=waveform_generator) # Run sampler. In this case we're going to use the `dynesty` sampler -result = bilby.core.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000, - injection_parameters=injection_parameters, outdir=outdir, label=label) +result = bilby.core.sampler.run_sampler( + likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000, + injection_parameters=injection_parameters, outdir=outdir, label=label) # make some plots of the outputs result.plot_corner() diff --git a/examples/injection_examples/using_gwin.py b/examples/injection_examples/using_gwin.py index a213c490c..dc312ddd8 100644 --- a/examples/injection_examples/using_gwin.py +++ b/examples/injection_examples/using_gwin.py @@ -48,7 +48,8 @@ fixed_parameters = injection_parameters.copy() for key in priors: fixed_parameters.pop(key) -# These lines generate the `model` object - see https://gwin.readthedocs.io/en/latest/api/gwin.models.gaussian_noise.html +# These lines generate the `model` object - see +# https://gwin.readthedocs.io/en/latest/api/gwin.models.gaussian_noise.html generator = FDomainDetFrameGenerator( FDomainCBCGenerator, 0., variable_args=variable_parameters, detectors=['H1', 'L1'], @@ -64,6 +65,7 @@ model.update(**injection_parameters) # This create a dummy class to convert the model into a bilby.likelihood object class GWINLikelihood(bilby.core.likelihood.Likelihood): + def __init__(self, model): """ A likelihood to wrap around GWIN model objects diff --git a/examples/other_examples/linear_regression_pymc3.py b/examples/other_examples/linear_regression_pymc3.py index c10d5d3c9..df98e5a5b 100644 --- a/examples/other_examples/linear_regression_pymc3.py +++ b/examples/other_examples/linear_regression_pymc3.py @@ -18,10 +18,12 @@ label = 'linear_regression_pymc3' outdir = 'outdir' bilby.utils.check_directory_exists_and_if_not_mkdir(outdir) + # First, we define our "signal model", in this case a simple linear function def model(time, m, c): return time * m + c + # Now we define the injection parameters which we make simulated data with injection_parameters = dict(m=0.5, c=0.2) @@ -51,7 +53,7 @@ likelihood = GaussianLikelihood(time, data, model, sigma=sigma) # From hereon, the syntax is exactly equivalent to other bilby examples # We make a prior -priors = {} +priors = dict() priors['m'] = bilby.core.prior.Uniform(0, 5, 'm') priors['c'] = bilby.core.prior.Uniform(-2, 2, 'c') diff --git a/examples/other_examples/linear_regression_pymc3_custom_likelihood.py b/examples/other_examples/linear_regression_pymc3_custom_likelihood.py index 9321ea4c8..f62b187ac 100644 --- a/examples/other_examples/linear_regression_pymc3_custom_likelihood.py +++ b/examples/other_examples/linear_regression_pymc3_custom_likelihood.py @@ -20,10 +20,12 @@ label = 'linear_regression_pymc3_custom_likelihood' outdir = 'outdir' bilby.utils.check_directory_exists_and_if_not_mkdir(outdir) + # First, we define our "signal model", in this case a simple linear function def model(time, m, c): return time * m + c + # Now we define the injection parameters which we make simulated data with injection_parameters = dict(m=0.5, c=0.2) @@ -51,6 +53,7 @@ fig.savefig('{}/{}_data.png'.format(outdir, label)) # Parameter estimation: we now define a Gaussian Likelihood class relevant for # our model. class GaussianLikelihoodPyMC3(bilby.Likelihood): + def __init__(self, x, y, sigma, function): """ A general Gaussian likelihood - the parameters are inferred from the @@ -105,6 +108,7 @@ class GaussianLikelihoodPyMC3(bilby.Likelihood): # set the likelihood distribution pm.Normal('likelihood', mu=mu, sd=self.sigma, observed=self.y) + # Now lets instantiate a version of our GaussianLikelihood, giving it # the time, data and signal model likelihood = GaussianLikelihoodPyMC3(time, data, sigma, model) @@ -144,6 +148,6 @@ priors['c'] = PriorPyMC3(-2, 2, 'c') # And run sampler result = bilby.run_sampler( likelihood=likelihood, priors=priors, sampler='pymc3', draws=1000, - tune=1000, discard_tuned_samples=True, injection_parameters=injection_parameters, - outdir=outdir, label=label) + tune=1000, discard_tuned_samples=True, + injection_parameters=injection_parameters, outdir=outdir, label=label) result.plot_corner() diff --git a/examples/supernova_example/supernova_example.py b/examples/supernova_example/supernova_example.py index b6eeb99a1..05b2eebcd 100644 --- a/examples/supernova_example/supernova_example.py +++ b/examples/supernova_example/supernova_example.py @@ -1,4 +1,4 @@ -#!/bin/python +#!/usr/bin/env python """ Tutorial to demonstrate running parameter estimation/model selection on an NR @@ -9,10 +9,10 @@ factor. (See https://arxiv.org/pdf/1202.3256.pdf) """ from __future__ import division, print_function import numpy as np +import bilby -# Set the duration and sampling frequency of the data segment that we're going to inject the signal into -import bilby.gw.likelihood - +# Set the duration and sampling frequency of the data segment that we're going +# to inject the signal into duration = 3. sampling_frequency = 4096. -- GitLab