diff --git a/CHANGELOG.md b/CHANGELOG.md index 1fdaae3edc4c09906e1af38f9702a232be97a9d1..a5f45836dce4a12637061f1cb292e722de210f1c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,7 +26,9 @@ Changes currently on master, but not under a tag. - Autocorrelation calculation moved into parent class - Fix interpretation of kwargs for dynesty - PowerSpectralDensity structure modified -- Fixed bug in get_open_data +- Fixed bug in get_open_data +- Modified how sampling in non-standard parameters is done, the + `non_standard_sampling_parameter_keys` kwarg has been removed - .prior files are no longer created. The prior is stored in the result object. - Changed the way repr works for priors. The repr can now be used to re-instantiate the Prior in most cases diff --git a/examples/injection_examples/basic_tutorial.py b/examples/injection_examples/basic_tutorial.py index a8575e17b28fb9cea657e979fb048edab8a00e65..2c18ef1c33fba11302dd09a6a892196a2e8861dd 100644 --- a/examples/injection_examples/basic_tutorial.py +++ b/examples/injection_examples/basic_tutorial.py @@ -40,18 +40,16 @@ waveform_generator = tupak.gw.WaveformGenerator( duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, parameters=injection_parameters, waveform_arguments=waveform_arguments) -hf_signal = waveform_generator.frequency_domain_strain() # Set up interferometers. In this case we'll use two interferometers # (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design # sensitivity -IFOs = [] -for name in ["H1", "L1"]: - IFOs.append( - tupak.gw.detector.get_interferometer_with_fake_noise_and_injection( - name, injection_polarizations=hf_signal, - injection_parameters=injection_parameters, duration=duration, - sampling_frequency=sampling_frequency, outdir=outdir)) +ifos = tupak.gw.detector.InterferometerList(['H1', 'L1']) +ifos.set_strain_data_from_power_spectral_densities( + sampling_frequency=sampling_frequency, duration=duration, + start_time=injection_parameters['geocent_time']-3) +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, @@ -68,7 +66,7 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi', 'ra', ' priors[key] = injection_parameters[key] # Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator -likelihood = tupak.gw.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator, +likelihood = tupak.gw.GravitationalWaveTransient(interferometers=ifos, waveform_generator=waveform_generator, time_marginalization=False, phase_marginalization=False, distance_marginalization=False, prior=priors) diff --git a/examples/injection_examples/change_sampled_parameters.py b/examples/injection_examples/change_sampled_parameters.py index 9f95dce802453bbd838f760db3d164b22a409b17..106ace05c782ec8dc485e46c288527722e8731f8 100644 --- a/examples/injection_examples/change_sampled_parameters.py +++ b/examples/injection_examples/change_sampled_parameters.py @@ -1,8 +1,10 @@ -#!/bin/python +#!/usr/bin/env python """ -Tutorial to demonstrate running parameter estimation sampling in non-standard parameters for an injected signal. +Tutorial to demonstrate running parameter estimation sampling in non-standard +parameters for an injected signal. -This example estimates the masses using a uniform prior in chirp mass, mass ratio and redshift. +This example estimates the masses using a uniform prior in chirp mass, +mass ratio and redshift. The cosmology is according to the Planck 2015 data release. """ from __future__ import division, print_function @@ -18,9 +20,10 @@ outdir = 'outdir' np.random.seed(151226) -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=500, iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413, - ra=1.375, dec=-1.2108) +injection_parameters = dict( + total_mass=66., mass_ratio=0.9, 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) waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', reference_frequency=50.) @@ -29,35 +32,45 @@ waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( sampling_frequency=sampling_frequency, duration=duration, frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, - parameter_conversion=tupak.gw.conversion.convert_to_lal_binary_black_hole_parameters, - non_standard_sampling_parameter_keys=['chirp_mass', 'mass_ratio', 'redshift'], - parameters=injection_parameters, waveform_arguments=waveform_arguments) -hf_signal = waveform_generator.frequency_domain_strain() + parameter_conversion= + tupak.gw.conversion.convert_to_lal_binary_black_hole_parameters, + waveform_arguments=waveform_arguments) # Set up interferometers. -IFOs = [tupak.gw.detector.get_interferometer_with_fake_noise_and_injection( - name, injection_polarizations=hf_signal, injection_parameters=injection_parameters, duration=duration, - sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']] +ifos = tupak.gw.detector.InterferometerList(['H1', 'L1', 'V1', 'K1']) +ifos.set_strain_data_from_power_spectral_densities( + sampling_frequency=sampling_frequency, duration=duration, + start_time=injection_parameters['geocent_time']-3) +ifos.inject_signal(waveform_generator=waveform_generator, + parameters=injection_parameters) # Set up prior priors = tupak.gw.prior.BBHPriorSet() 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) -priors['mass_ratio'] = tupak.prior.Uniform(name='mass_ratio', latex_label='$q$', minimum=0.1, maximum=1) -priors['redshift'] = tupak.prior.Uniform(name='redshift', latex_label='$z$', minimum=0, maximum=0.5) +priors['chirp_mass'] = tupak.prior.Uniform( + name='chirp_mass', latex_label='$m_c$', minimum=13, maximum=45) +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( + name='redshift', latex_label='$z$', minimum=0, maximum=0.5) # These parameters will not be sampled -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] +priors.pop('iota') +priors['cos_iota'] = np.cos(injection_parameters['iota']) print(priors) # Initialise GravitationalWaveTransient -likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator) +likelihood = tupak.gw.likelihood.GravitationalWaveTransient( + interferometers=ifos, waveform_generator=waveform_generator) # Run sampler -result = tupak.core.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', - injection_parameters=injection_parameters, label='DifferentParameters', - outdir=outdir, conversion_function=tupak.gw.conversion.generate_all_bbh_parameters) +result = tupak.core.sampler.run_sampler( + likelihood=likelihood, priors=priors, sampler='dynesty', outdir=outdir, + injection_parameters=injection_parameters, label='DifferentParameters', + conversion_function=tupak.gw.conversion.generate_all_bbh_parameters) result.plot_corner() diff --git a/examples/injection_examples/create_your_own_source_model.py b/examples/injection_examples/create_your_own_source_model.py index 604a513ef6bb9a70ba60a7db02e4d3dfecf6696f..8ee965ab9382370a71ea96ad42b22e882d7adc31 100644 --- a/examples/injection_examples/create_your_own_source_model.py +++ b/examples/injection_examples/create_your_own_source_model.py @@ -29,14 +29,14 @@ waveform_generator = tupak.gw.waveform_generator.WaveformGenerator(duration=dura sampling_frequency=sampling_frequency, frequency_domain_source_model=sine_gaussian, parameters=injection_parameters) -hf_signal = waveform_generator.frequency_domain_strain() # Set up interferometers. -IFOs = [tupak.gw.detector.get_interferometer_with_fake_noise_and_injection( - name, injection_polarizations=hf_signal, - injection_parameters=injection_parameters, duration=duration, - sampling_frequency=sampling_frequency, outdir=outdir) - for name in ['H1', 'L1', 'V1']] +ifos = tupak.gw.detector.InterferometerList(['H1', 'L1']) +ifos.set_strain_data_from_power_spectral_densities( + sampling_frequency=sampling_frequency, duration=duration, + start_time=injection_parameters['geocent_time'] - 3) +ifos.inject_signal(waveform_generator=waveform_generator, + parameters=injection_parameters) # Here we define the priors for the search. We use the injection parameters # except for the amplitude, f0, and geocent_time @@ -44,7 +44,7 @@ prior = injection_parameters.copy() prior['A'] = tupak.core.prior.PowerLaw(alpha=-1, minimum=1e-25, maximum=1e-21, name='A') prior['f0'] = tupak.core.prior.Uniform(90, 110, 'f') -likelihood = tupak.gw.likelihood.GravitationalWaveTransient(IFOs, waveform_generator) +likelihood = tupak.gw.likelihood.GravitationalWaveTransient(ifos, waveform_generator) result = tupak.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 be48be57ab9c23f8594da2751f03e5be586fb698..b1f814c6e9cf3bc7fe6f2c5146cdf2f10fed49d8 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 @@ -36,23 +36,13 @@ waveform = tupak.gw.waveform_generator.WaveformGenerator(duration=duration, samp time_domain_source_model=time_domain_damped_sinusoid, parameters=injection_parameters) -hf_signal = waveform.frequency_domain_strain() -#note we could plot the time domain signal with the following code -# import matplotlib.pyplot as plt -# plt.plot(waveform.time_array, waveform.time_domain_strain()['plus']) - -# or the frequency-domain signal: -# plt.loglog(waveform.frequency_array, abs(waveform.frequency_domain_strain()['plus'])) - - - # inject the signal into three interferometers -IFOs = [tupak.gw.detector.get_interferometer_with_fake_noise_and_injection( - name, injection_polarizations=hf_signal, - injection_parameters=injection_parameters, duration=duration, - sampling_frequency=sampling_frequency, outdir=outdir) - for name in ['H1', 'L1']] - +ifos = tupak.gw.detector.InterferometerList(['H1', 'L1']) +ifos.set_strain_data_from_power_spectral_densities( + sampling_frequency=sampling_frequency, duration=duration, + start_time=injection_parameters['geocent_time'] - 3) +ifos.inject_signal(waveform_generator=waveform, + parameters=injection_parameters) # create the priors prior = injection_parameters.copy() @@ -61,9 +51,8 @@ prior['damping_time'] = tupak.core.prior.Uniform(0, 1, r'damping time') prior['frequency'] = tupak.core.prior.Uniform(0, 200, r'frequency') prior['phase'] = tupak.core.prior.Uniform(-np.pi / 2, np.pi / 2, r'$\phi$') - # define likelihood -likelihood = tupak.gw.likelihood.GravitationalWaveTransient(IFOs, waveform) +likelihood = tupak.gw.likelihood.GravitationalWaveTransient(ifos, waveform) # launch sampler result = tupak.core.sampler.run_sampler(likelihood, prior, sampler='dynesty', npoints=1000, diff --git a/examples/injection_examples/eccentric_inspiral.py b/examples/injection_examples/eccentric_inspiral.py index 79bdf8088e8cee636b389132ff09a5e0d0ff1937..f4077b640c2d3446f7c32699e58e4bdfda3f3215 100644 --- a/examples/injection_examples/eccentric_inspiral.py +++ b/examples/injection_examples/eccentric_inspiral.py @@ -40,7 +40,6 @@ waveform_generator = tupak.gw.WaveformGenerator( frequency_domain_source_model=tupak.gw.source.lal_eccentric_binary_black_hole_no_spins, parameters=injection_parameters, waveform_arguments=waveform_arguments) -hf_signal = waveform_generator.frequency_domain_strain() # Setting up three interferometers (LIGO-Hanford (H1), LIGO-Livingston (L1), and # Virgo (V1)) at their design sensitivities. The maximum frequency is set just @@ -49,13 +48,15 @@ hf_signal = waveform_generator.frequency_domain_strain() minimum_frequency = 10. maximum_frequency = 128. -IFOs = tupak.gw.detector.InterferometerList(['H1', 'L1', 'V1']) -for IFO in IFOs: - IFO.minimum_frequency = minimum_frequency - IFO.maximum_frequency = maximum_frequency - -IFOs.set_strain_data_from_power_spectral_densities(sampling_frequency, duration) -IFOs.inject_signal(waveform_generator=waveform_generator, parameters=injection_parameters) +ifos = tupak.gw.detector.InterferometerList(['H1', 'L1']) +for ifo in ifos: + ifo.minimum_frequency = minimum_frequency + ifo.maximum_frequency = maximum_frequency +ifos.set_strain_data_from_power_spectral_densities( + sampling_frequency=sampling_frequency, duration=duration, + start_time=injection_parameters['geocent_time'] - 3) +ifos.inject_signal(waveform_generator=waveform_generator, + parameters=injection_parameters) # Now we set up the priors on each of the binary parameters. priors = dict() @@ -71,7 +72,7 @@ priors["phase"] = tupak.core.prior.Uniform(name='phase', minimum=0, maximum=2 * priors["geocent_time"] = tupak.core.prior.Uniform(1180002600.9, 1180002601.1, name='geocent_time') # Initialising the likelihood function. -likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers=IFOs, +likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers=ifos, waveform_generator=waveform_generator, time_marginalization=False, phase_marginalization=False, distance_marginalization=False, prior=priors) diff --git a/examples/injection_examples/how_to_specify_the_prior.py b/examples/injection_examples/how_to_specify_the_prior.py index 450250db63b3c141ca87e9ef1689328b278dd857..51b7877dde665ece58094500d25e9504e81c94d7 100644 --- a/examples/injection_examples/how_to_specify_the_prior.py +++ b/examples/injection_examples/how_to_specify_the_prior.py @@ -27,12 +27,14 @@ waveform_generator = tupak.gw.WaveformGenerator( duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, parameters=injection_parameters, waveform_arguments=waveform_arguments) -hf_signal = waveform_generator.frequency_domain_strain() # Set up interferometers. -IFOs = [tupak.gw.detector.get_interferometer_with_fake_noise_and_injection( - name, injection_polarizations=hf_signal, injection_parameters=injection_parameters, duration=duration, - sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']] +ifos = tupak.gw.detector.InterferometerList(['H1', 'L1']) +ifos.set_strain_data_from_power_spectral_densities( + sampling_frequency=sampling_frequency, duration=duration, + start_time=injection_parameters['geocent_time'] - 3) +ifos.inject_signal(waveform_generator=waveform_generator, + parameters=injection_parameters) # Set up prior # This loads in a predefined set of priors for BBHs. @@ -56,7 +58,7 @@ priors['a_2'] = tupak.core.prior.Interped(name='a_2', xx=a_2, yy=p_a_2, minimum= # Enjoy. # Initialise GravitationalWaveTransient -likelihood = tupak.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator) +likelihood = tupak.GravitationalWaveTransient(interferometers=ifos, waveform_generator=waveform_generator) # Run sampler result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', diff --git a/examples/injection_examples/marginalized_likelihood.py b/examples/injection_examples/marginalized_likelihood.py index b91cee3e109f257d55b53c8238b924dbef870d79..8c05d0f977af14b6cd19c9614c2be956db71dcf3 100644 --- a/examples/injection_examples/marginalized_likelihood.py +++ b/examples/injection_examples/marginalized_likelihood.py @@ -26,16 +26,14 @@ waveform_generator = tupak.gw.WaveformGenerator( duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, parameters=injection_parameters, waveform_arguments=waveform_arguments) -hf_signal = waveform_generator.frequency_domain_strain() # Set up interferometers. -IFOs = [] -for name in ["H1", "L1", "V1"]: - IFOs.append( - tupak.gw.detector.get_interferometer_with_fake_noise_and_injection( - name, injection_polarizations=hf_signal, - injection_parameters=injection_parameters, duration=duration, - sampling_frequency=sampling_frequency, outdir=outdir)) +ifos = tupak.gw.detector.InterferometerList(['H1', 'L1']) +ifos.set_strain_data_from_power_spectral_densities( + sampling_frequency=sampling_frequency, duration=duration, + start_time=injection_parameters['geocent_time'] - 3) +ifos.inject_signal(waveform_generator=waveform_generator, + parameters=injection_parameters) # Set up prior priors = tupak.gw.prior.BBHPriorSet() @@ -47,7 +45,7 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'iota', 'ra', # 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. likelihood = tupak.gw.GravitationalWaveTransient( - interferometers=IFOs, waveform_generator=waveform_generator, prior=priors, + interferometers=ifos, waveform_generator=waveform_generator, prior=priors, distance_marginalization=False, phase_marginalization=True, time_marginalization=False) diff --git a/examples/injection_examples/plot_time_domain_data.py b/examples/injection_examples/plot_time_domain_data.py index a0b1a61a889cafb1569030680fec478f12d4da5b..37a0868c84157b2261a57499b9985252fc63c4c5 100644 --- a/examples/injection_examples/plot_time_domain_data.py +++ b/examples/injection_examples/plot_time_domain_data.py @@ -25,7 +25,7 @@ waveform_generator = tupak.gw.WaveformGenerator( duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, parameters=injection_parameters, waveform_arguments=waveform_arguments) -hf_signal = waveform_generator.frequency_domain_strain() +hf_signal = waveform_generator.frequency_domain_strain(injection_parameters) H1 = tupak.gw.detector.get_interferometer_with_fake_noise_and_injection( 'H1', injection_polarizations=hf_signal, diff --git a/examples/injection_examples/sine_gaussian_example.py b/examples/injection_examples/sine_gaussian_example.py index 8761cf6b388542a34717ec52c1c71c637d6fac1f..aba4081847b49c4f7388df6ccf348c4558c8c5d7 100644 --- a/examples/injection_examples/sine_gaussian_example.py +++ b/examples/injection_examples/sine_gaussian_example.py @@ -29,13 +29,15 @@ waveform_generator = tupak.gw.waveform_generator.WaveformGenerator(duration=dura sampling_frequency=sampling_frequency, frequency_domain_source_model=tupak.gw.source.sinegaussian, parameters=injection_parameters) -hf_signal = waveform_generator.frequency_domain_strain() # 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 = [tupak.gw.detector.get_interferometer_with_fake_noise_and_injection( - name, injection_polarizations=hf_signal, injection_parameters=injection_parameters, duration=duration, - sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']] +ifos = tupak.gw.detector.InterferometerList(['H1', 'L1', 'V1']) +ifos.set_strain_data_from_power_spectral_densities( + sampling_frequency=sampling_frequency, duration=duration, + start_time=injection_parameters['geocent_time']-3) +ifos.inject_signal(waveform_generator=waveform_generator, + parameters=injection_parameters) # Set up prior, which is a dictionary priors = dict() @@ -55,7 +57,7 @@ priors['frequency'] = tupak.core.prior.Uniform(30, 1000, 'frequency') 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 -likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator) +likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers=ifos, waveform_generator=waveform_generator) # Run sampler. In this case we're going to use the `dynesty` sampler result = tupak.core.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000, diff --git a/examples/supernova_example/supernova_example.py b/examples/supernova_example/supernova_example.py index 86ba46c80afb207a4833570666c9eb14417edf7f..2365bd55ef652f47819214acf7635615e5ab82b2 100644 --- a/examples/supernova_example/supernova_example.py +++ b/examples/supernova_example/supernova_example.py @@ -13,7 +13,7 @@ import numpy as np # Set the duration and sampling frequency of the data segment that we're going to inject the signal into import tupak.gw.likelihood -time_duration = 3. +duration = 3. sampling_frequency = 4096. # Specify the output directory and the name of the simulation. @@ -35,19 +35,19 @@ injection_parameters = dict(file_path='MuellerL15_example_inj.txt', # Create the waveform_generator using a supernova source function waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( - time_duration=time_duration, sampling_frequency=sampling_frequency, + duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=tupak.gw.source.supernova, parameters=injection_parameters) -hf_signal = waveform_generator.frequency_domain_strain() # 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 = [tupak.gw.detector.get_interferometer_with_fake_noise_and_injection( - name, injection_polarizations=hf_signal, - injection_parameters=injection_parameters, time_duration=time_duration, - sampling_frequency=sampling_frequency, outdir=outdir) - for name in ['H1', 'L1', 'V1']] +ifos = tupak.gw.detector.InterferometerList(['H1', 'L1']) +ifos.set_strain_data_from_power_spectral_densities( + sampling_frequency=sampling_frequency, duration=duration, + start_time=injection_parameters['geocent_time']-3) +ifos.inject_signal(waveform_generator=waveform_generator, + parameters=injection_parameters) # read in from a file the PCs used to create the signal model. realPCs = np.loadtxt('SupernovaRealPCs.txt') @@ -59,7 +59,7 @@ simulation_parameters = dict( realPCs=realPCs, imagPCs=imagPCs) search_waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( - time_duration=time_duration, sampling_frequency=sampling_frequency, + duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=tupak.gw.source.supernova_pca_model, waveform_arguments=simulation_parameters) @@ -85,7 +85,7 @@ priors['geocent_time'] = tupak.core.prior.Uniform( # Initialise the likelihood by passing in the interferometer data (IFOs) and # the waveoform generator likelihood = tupak.gw.likelihood.GravitationalWaveTransient( - interferometers=IFOs, waveform_generator=search_waveform_generator) + interferometers=ifos, waveform_generator=search_waveform_generator) # Run sampler. result = tupak.run_sampler( diff --git a/test/conversion_test.py b/test/conversion_test.py index f0e4409f2d360c93ba93a96f6ac33c7b4089dc88..6a57a177359fd08c92df411ebf9551745c51513f 100644 --- a/test/conversion_test.py +++ b/test/conversion_test.py @@ -1,7 +1,7 @@ from __future__ import division, absolute_import import unittest import mock -import tupak +from tupak.gw import conversion import numpy as np @@ -47,52 +47,52 @@ class TestBasicConversions(unittest.TestCase): del self.symmetric_mass_ratio def test_total_mass_and_mass_ratio_to_component_masses(self): - mass_1, mass_2 = tupak.gw.conversion.total_mass_and_mass_ratio_to_component_masses(self.mass_ratio, self.total_mass) + mass_1, mass_2 = conversion.total_mass_and_mass_ratio_to_component_masses(self.mass_ratio, self.total_mass) self.assertTrue(all([abs(mass_1 - self.mass_1) < 1e-5, abs(mass_2 - self.mass_2) < 1e-5])) def test_symmetric_mass_ratio_to_mass_ratio(self): - mass_ratio = tupak.gw.conversion.symmetric_mass_ratio_to_mass_ratio(self.symmetric_mass_ratio) + mass_ratio = conversion.symmetric_mass_ratio_to_mass_ratio(self.symmetric_mass_ratio) self.assertAlmostEqual(self.mass_ratio, mass_ratio) def test_chirp_mass_and_total_mass_to_symmetric_mass_ratio(self): - symmetric_mass_ratio = tupak.gw.conversion.chirp_mass_and_total_mass_to_symmetric_mass_ratio(self.chirp_mass, self.total_mass) + symmetric_mass_ratio = conversion.chirp_mass_and_total_mass_to_symmetric_mass_ratio(self.chirp_mass, self.total_mass) self.assertAlmostEqual(self.symmetric_mass_ratio, symmetric_mass_ratio) def test_chirp_mass_and_mass_ratio_to_total_mass(self): - total_mass = tupak.gw.conversion.chirp_mass_and_mass_ratio_to_total_mass(self.chirp_mass, self.mass_ratio) + total_mass = conversion.chirp_mass_and_mass_ratio_to_total_mass(self.chirp_mass, self.mass_ratio) self.assertAlmostEqual(self.total_mass, total_mass) def test_component_masses_to_chirp_mass(self): - chirp_mass = tupak.gw.conversion.component_masses_to_chirp_mass(self.mass_1, self.mass_2) + chirp_mass = conversion.component_masses_to_chirp_mass(self.mass_1, self.mass_2) self.assertAlmostEqual(self.chirp_mass, chirp_mass) def test_component_masses_to_total_mass(self): - total_mass = tupak.gw.conversion.component_masses_to_total_mass(self.mass_1, self.mass_2) + total_mass = conversion.component_masses_to_total_mass(self.mass_1, self.mass_2) self.assertAlmostEqual(self.total_mass, total_mass) def test_component_masses_to_symmetric_mass_ratio(self): - symmetric_mass_ratio = tupak.gw.conversion.component_masses_to_symmetric_mass_ratio(self.mass_1, self.mass_2) + symmetric_mass_ratio = conversion.component_masses_to_symmetric_mass_ratio(self.mass_1, self.mass_2) self.assertAlmostEqual(self.symmetric_mass_ratio, symmetric_mass_ratio) def test_component_masses_to_mass_ratio(self): - mass_ratio = tupak.gw.conversion.component_masses_to_mass_ratio(self.mass_1, self.mass_2) + mass_ratio = conversion.component_masses_to_mass_ratio(self.mass_1, self.mass_2) self.assertAlmostEqual(self.mass_ratio, mass_ratio) def test_mass_1_and_chirp_mass_to_mass_ratio(self): - mass_ratio = tupak.gw.conversion.mass_1_and_chirp_mass_to_mass_ratio(self.mass_1, self.chirp_mass) + mass_ratio = conversion.mass_1_and_chirp_mass_to_mass_ratio(self.mass_1, self.chirp_mass) self.assertAlmostEqual(self.mass_ratio, mass_ratio) def test_lambda_tilde_to_lambda_1_lambda_2(self): lambda_1, lambda_2 =\ - tupak.gw.conversion.lambda_tilde_to_lambda_1_lambda_2( + conversion.lambda_tilde_to_lambda_1_lambda_2( self.lambda_tilde, self.mass_1, self.mass_2) self.assertTrue(all([abs(self.lambda_1 - lambda_1) < 1e-5, abs(self.lambda_2 - lambda_2) < 1e-5])) def test_lambda_tilde_delta_lambda_to_lambda_1_lambda_2(self): lambda_1, lambda_2 =\ - tupak.gw.conversion.lambda_tilde_delta_lambda_to_lambda_1_lambda_2( + conversion.lambda_tilde_delta_lambda_to_lambda_1_lambda_2( self.lambda_tilde, self.delta_lambda, self.mass_1, self.mass_2) self.assertTrue(all([abs(self.lambda_1 - lambda_1) < 1e-5, abs(self.lambda_2 - lambda_2) < 1e-5])) @@ -113,17 +113,10 @@ class TestConvertToLALBBHParams(unittest.TestCase): def test_cos_angle_to_angle_conversion(self): with mock.patch('numpy.arccos') as m: m.return_value = 42 - self.search_keys.append('cos_tilt_1') self.parameters['cos_tilt_1'] = 1 - self.parameters, _ = tupak.gw.conversion.convert_to_lal_binary_black_hole_parameters(self.parameters, self.search_keys) - self.assertEqual(42, self.parameters['tilt_1']) - - def test_cos_angle_to_angle_conversion_removal(self): - with mock.patch('numpy.arccos') as m: - m.return_value = 42 - self.search_keys.append('cos_tilt_1') - self.parameters['cos_tilt_1'] = 1 - self.parameters, _ = tupak.gw.conversion.convert_to_lal_binary_black_hole_parameters(self.parameters, self.search_keys, remove=True) + self.parameters, _ =\ + conversion.convert_to_lal_binary_black_hole_parameters( + self.parameters) self.assertDictEqual(self.parameters, dict(tilt_1=42, cos_tilt_1=1)) diff --git a/test/detector_test.py b/test/detector_test.py index a69ca162378220c3149148d31784c46972f2f69c..e9ff0e11af6dfa47f480c1e412574924dcb186b4 100644 --- a/test/detector_test.py +++ b/test/detector_test.py @@ -723,16 +723,6 @@ class TestInterferometerList(unittest.TestCase): self.ifo_list.inject_signal(parameters=None, waveform_generator=waveform_generator) self.assertTrue(m.called) - @patch.object(tupak.gw.waveform_generator.WaveformGenerator, 'frequency_domain_strain') - def test_inject_signal_pol_none_sets_wg_parameters(self, m): - waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( - frequency_domain_source_model=lambda x, y, z: x) - parameters = dict(y=1, z=2) - self.ifo1.inject_signal = MagicMock(return_value=None) - self.ifo2.inject_signal = MagicMock(return_value=None) - self.ifo_list.inject_signal(parameters=parameters, waveform_generator=waveform_generator) - self.assertDictEqual(parameters, waveform_generator.parameters) - @patch.object(tupak.gw.detector.Interferometer, 'inject_signal') def test_inject_signal_with_inj_pol(self, m): self.ifo_list.inject_signal(injection_polarizations=dict(plus=1)) diff --git a/test/gw_likelihood_test.py b/test/gw_likelihood_test.py index fa8e53b6ec255b5defd1ad98ccd7da0f9f39b46b..e04f16a453291684ceb620d6e5346bbd210a6ec4 100644 --- a/test/gw_likelihood_test.py +++ b/test/gw_likelihood_test.py @@ -1,9 +1,7 @@ from __future__ import division, absolute_import import unittest -import mock import tupak import numpy as np -from scipy.special import logsumexp class TestBasicGWTransient(unittest.TestCase): @@ -19,7 +17,7 @@ class TestBasicGWTransient(unittest.TestCase): self.interferometers.set_strain_data_from_power_spectral_densities( sampling_frequency=2048, duration=4) self.waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( - duration=4, sampling_frequency=2048, parameters=self.parameters, + duration=4, sampling_frequency=2048, frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, ) @@ -27,6 +25,7 @@ class TestBasicGWTransient(unittest.TestCase): interferometers=self.interferometers, waveform_generator=self.waveform_generator ) + self.likelihood.parameters = self.parameters.copy() def tearDown(self): del self.parameters @@ -56,10 +55,10 @@ class TestBasicGWTransient(unittest.TestCase): def test_likelihood_zero_when_waveform_is_none(self): """Test log likelihood returns np.nan_to_num(-np.inf) when the waveform is None""" - self.likelihood.waveform_generator.parameters['mass_2'] = 32 + self.likelihood.parameters['mass_2'] = 32 self.assertEqual(self.likelihood.log_likelihood_ratio(), np.nan_to_num(-np.inf)) - self.likelihood.waveform_generator.parameters['mass_2'] = 29 + self.likelihood.parameters['mass_2'] = 29 def test_repr(self): expected = 'BasicGravitationalWaveTransient(interferometers={},\n\twaveform_generator={})'.format( @@ -83,7 +82,6 @@ class TestGWTransient(unittest.TestCase): sampling_frequency=self.sampling_frequency, duration=self.duration) self.waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( duration=self.duration, sampling_frequency=self.sampling_frequency, - parameters=self.parameters.copy(), frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, ) @@ -96,12 +94,7 @@ class TestGWTransient(unittest.TestCase): interferometers=self.interferometers, waveform_generator=self.waveform_generator, prior=self.prior.copy() ) - - # self.distance = tupak.gw.likelihood.GravitationalWaveTransient( - # interferometers=self.interferometers, - # waveform_generator=self.waveform_generator, - # distance_marginalization=True, prior=self.prior - # ) + self.likelihood.parameters = self.parameters.copy() def tearDown(self): del self.parameters @@ -109,7 +102,6 @@ class TestGWTransient(unittest.TestCase): del self.waveform_generator del self.prior del self.likelihood - # del self.distance def test_noise_log_likelihood(self): """Test noise log likelihood matches precomputed value""" @@ -133,10 +125,10 @@ class TestGWTransient(unittest.TestCase): def test_likelihood_zero_when_waveform_is_none(self): """Test log likelihood returns np.nan_to_num(-np.inf) when the waveform is None""" - self.likelihood.waveform_generator.parameters['mass_2'] = 32 + self.likelihood.parameters['mass_2'] = 32 self.assertEqual(self.likelihood.log_likelihood_ratio(), np.nan_to_num(-np.inf)) - self.likelihood.waveform_generator.parameters['mass_2'] = 29 + self.likelihood.parameters['mass_2'] = 29 def test_repr(self): expected = 'GravitationalWaveTransient(interferometers={},\n\twaveform_generator={},\n\t' \ @@ -163,7 +155,6 @@ class TestTimeMarginalization(unittest.TestCase): self.waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( duration=self.duration, sampling_frequency=self.sampling_frequency, - parameters=self.parameters.copy(), frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, ) @@ -182,6 +173,8 @@ class TestTimeMarginalization(unittest.TestCase): waveform_generator=self.waveform_generator, time_marginalization=True, prior=self.prior.copy() ) + for like in [self.likelihood, self.time]: + like.parameters = self.parameters.copy() def tearDown(self): del self.duration @@ -199,12 +192,12 @@ class TestTimeMarginalization(unittest.TestCase): times = np.linspace(self.prior['geocent_time'].minimum, self.prior['geocent_time'].maximum, 4097)[:-1] for time in times: - self.waveform_generator.parameters['geocent_time'] = time + self.likelihood.parameters['geocent_time'] = time like.append(np.exp(self.likelihood.log_likelihood_ratio())) marg_like = np.log(np.trapz(like, times) / self.waveform_generator.duration) - self.waveform_generator.parameters = self.parameters.copy() + self.time.parameters = self.parameters.copy() self.assertAlmostEqual(marg_like, self.time.log_likelihood_ratio(), delta=0.5) @@ -228,7 +221,6 @@ class TestMarginalizedLikelihood(unittest.TestCase): self.waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( duration=self.duration, sampling_frequency=self.sampling_frequency, - parameters=self.parameters.copy(), frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, ) @@ -292,7 +284,6 @@ class TestPhaseMarginalization(unittest.TestCase): self.waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( duration=self.duration, sampling_frequency=self.sampling_frequency, - parameters=self.parameters.copy(), frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, ) @@ -311,6 +302,8 @@ class TestPhaseMarginalization(unittest.TestCase): waveform_generator=self.waveform_generator, phase_marginalization=True, prior=self.prior.copy() ) + for like in [self.likelihood, self.phase]: + like.parameters = self.parameters.copy() def tearDown(self): del self.duration @@ -327,11 +320,11 @@ class TestPhaseMarginalization(unittest.TestCase): like = [] phases = np.linspace(0, 2 * np.pi, 1000) for phase in phases: - self.waveform_generator.parameters['phase'] = phase + self.likelihood.parameters['phase'] = phase like.append(np.exp(self.likelihood.log_likelihood_ratio())) marg_like = np.log(np.trapz(like, phases) / (2 * np.pi)) - self.waveform_generator.parameters = self.parameters.copy() + self.phase.parameters = self.parameters.copy() self.assertAlmostEqual(marg_like, self.phase.log_likelihood_ratio(), delta=0.5) @@ -354,7 +347,6 @@ class TestTimePhaseMarginalization(unittest.TestCase): self.waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( duration=self.duration, sampling_frequency=self.sampling_frequency, - parameters=self.parameters.copy(), frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, ) @@ -386,6 +378,8 @@ class TestTimePhaseMarginalization(unittest.TestCase): time_marginalization=True, phase_marginalization=True, prior=self.prior.copy() ) + for like in [self.likelihood, self.time, self.phase, self.time_phase]: + like.parameters = self.parameters.copy() def tearDown(self): del self.duration @@ -405,12 +399,12 @@ class TestTimePhaseMarginalization(unittest.TestCase): times = np.linspace(self.prior['geocent_time'].minimum, self.prior['geocent_time'].maximum, 4097)[:-1] for time in times: - self.waveform_generator.parameters['geocent_time'] = time + self.phase.parameters['geocent_time'] = time like.append(np.exp(self.phase.log_likelihood_ratio())) marg_like = np.log(np.trapz(like, times) / self.waveform_generator.duration) - self.waveform_generator.parameters = self.parameters.copy() + self.time_phase.parameters = self.parameters.copy() self.assertAlmostEqual(marg_like, self.time_phase.log_likelihood_ratio(), delta=0.5) @@ -418,11 +412,11 @@ class TestTimePhaseMarginalization(unittest.TestCase): like = [] phases = np.linspace(0, 2 * np.pi, 1000) for phase in phases: - self.waveform_generator.parameters['phase'] = phase + self.time.parameters['phase'] = phase like.append(np.exp(self.time.log_likelihood_ratio())) marg_like = np.log(np.trapz(like, phases) / (2 * np.pi)) - self.waveform_generator.parameters = self.parameters.copy() + self.time_phase.parameters = self.parameters.copy() self.assertAlmostEqual(marg_like, self.time_phase.log_likelihood_ratio(), delta=0.5) @@ -435,7 +429,6 @@ class TestBBHLikelihoodSetUp(unittest.TestCase): def tearDown(self): del self.ifos - del self.like def test_instantiation(self): self.like = tupak.gw.likelihood.get_binary_black_hole_likelihood( diff --git a/test/waveform_generator_test.py b/test/waveform_generator_test.py index a7cd95ebe9ee6ce43f4bbe5dbbf5a04bf44f9c2a..a03f250a3db219bbc810c8da700e007ad65ceaa1 100644 --- a/test/waveform_generator_test.py +++ b/test/waveform_generator_test.py @@ -20,8 +20,9 @@ class TestWaveformGeneratorInstantiationWithoutOptionalParameters(unittest.TestC def setUp(self): self.waveform_generator = \ - tupak.gw.waveform_generator.WaveformGenerator(1, 4096, - frequency_domain_source_model=dummy_func_dict_return_value) + tupak.gw.waveform_generator.WaveformGenerator( + 1, 4096, + frequency_domain_source_model=dummy_func_dict_return_value) self.simulation_parameters = dict(amplitude=1e-21, mu=100, sigma=1, ra=1.375, dec=-1.2108, @@ -34,16 +35,14 @@ class TestWaveformGeneratorInstantiationWithoutOptionalParameters(unittest.TestC def test_repr(self): expected = 'WaveformGenerator(duration={}, sampling_frequency={}, start_time={}, ' \ - 'frequency_domain_source_model={}, time_domain_source_model={}, parameters={}, ' \ - 'parameter_conversion={}, non_standard_sampling_parameter_keys={}, waveform_arguments={})'\ + 'frequency_domain_source_model={}, time_domain_source_model={}, ' \ + 'parameter_conversion={}, waveform_arguments={})'\ .format(self.waveform_generator.duration, self.waveform_generator.sampling_frequency, self.waveform_generator.start_time, self.waveform_generator.frequency_domain_source_model.__name__, self.waveform_generator.time_domain_source_model, - self.waveform_generator.parameters, None, - self.waveform_generator.non_standard_sampling_parameter_keys, self.waveform_generator.waveform_arguments) self.assertEqual(expected, repr(self.waveform_generator)) @@ -52,16 +51,14 @@ class TestWaveformGeneratorInstantiationWithoutOptionalParameters(unittest.TestC tupak.gw.waveform_generator.WaveformGenerator(1, 4096, time_domain_source_model=dummy_func_dict_return_value) expected = 'WaveformGenerator(duration={}, sampling_frequency={}, start_time={}, ' \ - 'frequency_domain_source_model={}, time_domain_source_model={}, parameters={}, ' \ - 'parameter_conversion={}, non_standard_sampling_parameter_keys={}, waveform_arguments={})'\ + 'frequency_domain_source_model={}, time_domain_source_model={}, ' \ + 'parameter_conversion={}, waveform_arguments={})'\ .format(self.waveform_generator.duration, self.waveform_generator.sampling_frequency, self.waveform_generator.start_time, self.waveform_generator.frequency_domain_source_model, self.waveform_generator.time_domain_source_model.__name__, - self.waveform_generator.parameters, None, - self.waveform_generator.non_standard_sampling_parameter_keys, self.waveform_generator.waveform_arguments) self.assertEqual(expected, repr(self.waveform_generator)) @@ -71,16 +68,14 @@ class TestWaveformGeneratorInstantiationWithoutOptionalParameters(unittest.TestC self.waveform_generator.parameter_conversion = conversion_func expected = 'WaveformGenerator(duration={}, sampling_frequency={}, start_time={}, ' \ - 'frequency_domain_source_model={}, time_domain_source_model={}, parameters={}, ' \ - 'parameter_conversion={}, non_standard_sampling_parameter_keys={}, waveform_arguments={})'\ + 'frequency_domain_source_model={}, time_domain_source_model={}, ' \ + 'parameter_conversion={}, waveform_arguments={})'\ .format(self.waveform_generator.duration, self.waveform_generator.sampling_frequency, self.waveform_generator.start_time, self.waveform_generator.frequency_domain_source_model.__name__, self.waveform_generator.time_domain_source_model, - self.waveform_generator.parameters, conversion_func.__name__, - self.waveform_generator.non_standard_sampling_parameter_keys, self.waveform_generator.waveform_arguments) self.assertEqual(expected, repr(self.waveform_generator)) @@ -100,6 +95,7 @@ class TestWaveformGeneratorInstantiationWithoutOptionalParameters(unittest.TestC self.assertIsInstance(self.waveform_generator.time_array, np.ndarray) def test_source_model_parameters(self): + self.waveform_generator.parameters = self.simulation_parameters.copy() self.assertListEqual(sorted(list(self.waveform_generator.parameters.keys())), sorted(list(self.simulation_parameters.keys()))) @@ -136,14 +132,15 @@ class TestSetters(unittest.TestCase): del self.simulation_parameters def test_parameter_setter_sets_expected_values_with_expected_keys(self): - self.waveform_generator.parameters = self.simulation_parameters + self.waveform_generator.parameters = self.simulation_parameters.copy() for key in self.simulation_parameters: self.assertEqual(self.waveform_generator.parameters[key], self.simulation_parameters[key]) def test_parameter_setter_none_handling(self): - self.waveform_generator.parameters = None - self.assertListEqual(sorted(list(self.waveform_generator.parameters.keys())), - sorted(list(self.simulation_parameters.keys()))) + with self.assertRaises(TypeError): + self.waveform_generator.parameters = None + # self.assertListEqual(sorted(list(self.waveform_generator.parameters.keys())), + # sorted(list(self.simulation_parameters.keys()))) def test_frequency_array_setter(self): new_frequency_array = np.arange(1, 100) @@ -157,11 +154,13 @@ class TestSetters(unittest.TestCase): def test_parameters_set_from_frequency_domain_source_model(self): self.waveform_generator.frequency_domain_source_model = dummy_func_dict_return_value + self.waveform_generator.parameters = self.simulation_parameters.copy() self.assertListEqual(sorted(list(self.waveform_generator.parameters.keys())), sorted(list(self.simulation_parameters.keys()))) def test_parameters_set_from_time_domain_source_model(self): self.waveform_generator.time_domain_source_model = dummy_func_dict_return_value + self.waveform_generator.parameters = self.simulation_parameters.copy() self.assertListEqual(sorted(list(self.waveform_generator.parameters.keys())), sorted(list(self.simulation_parameters.keys()))) @@ -195,10 +194,10 @@ class TestFrequencyDomainStrainMethod(unittest.TestCase): def test_parameter_conversion_is_called(self): self.waveform_generator.parameter_conversion = MagicMock(side_effect=KeyError('test')) with self.assertRaises(KeyError): - self.waveform_generator.frequency_domain_strain() + self.waveform_generator.frequency_domain_strain( + parameters=self.simulation_parameters) def test_frequency_domain_source_model_call(self): - self.waveform_generator.parameters = self.simulation_parameters expected = self.waveform_generator.frequency_domain_source_model(self.waveform_generator.frequency_array, self.simulation_parameters['amplitude'], self.simulation_parameters['mu'], @@ -207,36 +206,39 @@ class TestFrequencyDomainStrainMethod(unittest.TestCase): self.simulation_parameters['dec'], self.simulation_parameters['geocent_time'], self.simulation_parameters['psi']) - actual = self.waveform_generator.frequency_domain_strain() + actual = self.waveform_generator.frequency_domain_strain( + parameters=self.simulation_parameters) self.assertTrue(np.array_equal(expected['plus'], actual['plus'])) self.assertTrue(np.array_equal(expected['cross'], actual['cross'])) def test_time_domain_source_model_call_with_ndarray(self): self.waveform_generator.frequency_domain_source_model = None self.waveform_generator.time_domain_source_model = dummy_func_array_return_value - self.waveform_generator.parameters = self.simulation_parameters def side_effect(value, value2): return value with mock.patch('tupak.core.utils.nfft') as m: m.side_effect = side_effect - expected = self.waveform_generator.time_domain_strain() - actual = self.waveform_generator.frequency_domain_strain() + expected = self.waveform_generator.time_domain_strain( + parameters=self.simulation_parameters) + actual = self.waveform_generator.frequency_domain_strain( + parameters=self.simulation_parameters) self.assertTrue(np.array_equal(expected, actual)) def test_time_domain_source_model_call_with_dict(self): self.waveform_generator.frequency_domain_source_model = None self.waveform_generator.time_domain_source_model = dummy_func_dict_return_value - self.waveform_generator.parameters = self.simulation_parameters def side_effect(value, value2): return value, self.waveform_generator.frequency_array with mock.patch('tupak.core.utils.nfft') as m: m.side_effect = side_effect - expected = self.waveform_generator.time_domain_strain() - actual = self.waveform_generator.frequency_domain_strain() + expected = self.waveform_generator.time_domain_strain( + parameters=self.simulation_parameters) + actual = self.waveform_generator.frequency_domain_strain( + parameters=self.simulation_parameters) self.assertTrue(np.array_equal(expected['plus'], actual['plus'])) self.assertTrue(np.array_equal(expected['cross'], actual['cross'])) @@ -244,7 +246,8 @@ class TestFrequencyDomainStrainMethod(unittest.TestCase): self.waveform_generator.time_domain_source_model = None self.waveform_generator.frequency_domain_source_model = None with self.assertRaises(RuntimeError): - self.waveform_generator.frequency_domain_strain() + self.waveform_generator.frequency_domain_strain( + parameters=self.simulation_parameters) def test_key_popping(self): self.waveform_generator.parameter_conversion = MagicMock(return_value=(dict(amplitude=1e-21, mu=100, sigma=1, @@ -253,7 +256,8 @@ class TestFrequencyDomainStrainMethod(unittest.TestCase): psi=2.659, c=None, d=None), ['c', 'd'])) try: - self.waveform_generator.frequency_domain_strain() + self.waveform_generator.frequency_domain_strain( + parameters=self.simulation_parameters) except RuntimeError: pass self.assertListEqual(sorted(self.waveform_generator.parameters.keys()), @@ -279,10 +283,10 @@ class TestTimeDomainStrainMethod(unittest.TestCase): def test_parameter_conversion_is_called(self): self.waveform_generator.parameter_conversion = MagicMock(side_effect=KeyError('test')) with self.assertRaises(KeyError): - self.waveform_generator.time_domain_strain() + self.waveform_generator.time_domain_strain( + parameters=self.simulation_parameters) def test_time_domain_source_model_call(self): - self.waveform_generator.parameters = self.simulation_parameters expected = self.waveform_generator.time_domain_source_model(self.waveform_generator.time_array, self.simulation_parameters['amplitude'], self.simulation_parameters['mu'], @@ -291,36 +295,39 @@ class TestTimeDomainStrainMethod(unittest.TestCase): self.simulation_parameters['dec'], self.simulation_parameters['geocent_time'], self.simulation_parameters['psi']) - actual = self.waveform_generator.time_domain_strain() + actual = self.waveform_generator.time_domain_strain( + parameters=self.simulation_parameters) self.assertTrue(np.array_equal(expected['plus'], actual['plus'])) self.assertTrue(np.array_equal(expected['cross'], actual['cross'])) def test_frequency_domain_source_model_call_with_ndarray(self): self.waveform_generator.time_domain_source_model = None self.waveform_generator.frequency_domain_source_model = dummy_func_array_return_value - self.waveform_generator.parameters = self.simulation_parameters def side_effect(value, value2): return value with mock.patch('tupak.core.utils.infft') as m: m.side_effect = side_effect - expected = self.waveform_generator.frequency_domain_strain() - actual = self.waveform_generator.time_domain_strain() + expected = self.waveform_generator.frequency_domain_strain( + parameters=self.simulation_parameters) + actual = self.waveform_generator.time_domain_strain( + parameters=self.simulation_parameters) self.assertTrue(np.array_equal(expected, actual)) def test_frequency_domain_source_model_call_with_dict(self): self.waveform_generator.time_domain_source_model = None self.waveform_generator.frequency_domain_source_model = dummy_func_dict_return_value - self.waveform_generator.parameters = self.simulation_parameters def side_effect(value, value2): return value with mock.patch('tupak.core.utils.infft') as m: m.side_effect = side_effect - expected = self.waveform_generator.frequency_domain_strain() - actual = self.waveform_generator.time_domain_strain() + expected = self.waveform_generator.frequency_domain_strain( + parameters=self.simulation_parameters) + actual = self.waveform_generator.time_domain_strain( + parameters=self.simulation_parameters) self.assertTrue(np.array_equal(expected['plus'], actual['plus'])) self.assertTrue(np.array_equal(expected['cross'], actual['cross'])) @@ -328,7 +335,8 @@ class TestTimeDomainStrainMethod(unittest.TestCase): self.waveform_generator.time_domain_source_model = None self.waveform_generator.frequency_domain_source_model = None with self.assertRaises(RuntimeError): - self.waveform_generator.time_domain_strain() + self.waveform_generator.time_domain_strain( + parameters=self.simulation_parameters) def test_key_popping(self): self.waveform_generator.parameter_conversion = MagicMock(return_value=(dict(amplitude=1e-2, @@ -339,7 +347,8 @@ class TestTimeDomainStrainMethod(unittest.TestCase): psi=2.659, c=None, d=None), ['c', 'd'])) try: - self.waveform_generator.time_domain_strain() + self.waveform_generator.time_domain_strain( + parameters=self.simulation_parameters) except RuntimeError: pass self.assertListEqual(sorted(self.waveform_generator.parameters.keys()), diff --git a/tupak/core/prior.py b/tupak/core/prior.py index cfa01fb2e78f287a72643f5aa32ff79dc69afd39..3cc6f12f6b85150e52a358444c91477893b25ca7 100644 --- a/tupak/core/prior.py +++ b/tupak/core/prior.py @@ -135,12 +135,6 @@ class PriorSet(OrderedDict): 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: - if parameter in self: - continue - self[parameter] = create_default_prior(parameter, default_priors_file) - for missing_key in missing_keys: if not self.test_redundancy(missing_key): default_prior = create_default_prior(missing_key, default_priors_file) diff --git a/tupak/gw/conversion.py b/tupak/gw/conversion.py index 28c5b057da7fc99ad74d53729fb53db4c4e41f06..5f78db028e631ec5d6e0daf2ac8ffb221a48c9aa 100644 --- a/tupak/gw/conversion.py +++ b/tupak/gw/conversion.py @@ -47,7 +47,7 @@ def luminosity_distance_to_comoving_distance(distance): return redshift_to_comoving_distance(redshift) -def convert_to_lal_binary_black_hole_parameters(parameters, search_keys, remove=True): +def convert_to_lal_binary_black_hole_parameters(parameters): """ Convert parameters we have into parameters we need. @@ -65,10 +65,6 @@ def convert_to_lal_binary_black_hole_parameters(parameters, search_keys, remove= ---------- parameters: dict dictionary of parameter values to convert into the required parameters - search_keys: list - parameters which are needed for the waveform generation - remove: bool, optional - Whether or not to remove the extra key, necessary for sampling, default=True. Return ------ @@ -78,100 +74,89 @@ def convert_to_lal_binary_black_hole_parameters(parameters, search_keys, remove= keys which are added to parameters during function call """ - added_keys = [] converted_parameters = parameters.copy() - - if 'mass_1' not in search_keys and 'mass_2' not in search_keys: - if 'chirp_mass' in converted_parameters.keys(): - if 'total_mass' in converted_parameters.keys() and 'total_mass' not in added_keys: - converted_parameters['symmetric_mass_ratio'] = chirp_mass_and_total_mass_to_symmetric_mass_ratio( - converted_parameters['chirp_mass'], converted_parameters['total_mass']) - if remove: - added_keys.append('chirp_mass') - if 'symmetric_mass_ratio' in converted_parameters.keys() and 'symmetric_mass_ratio' not in added_keys: - converted_parameters['mass_ratio'] = \ - symmetric_mass_ratio_to_mass_ratio(converted_parameters['symmetric_mass_ratio']) - if remove: - added_keys.append('symmetric_mass_ratio') - if 'mass_ratio' in converted_parameters.keys() and 'mass_ratio' not in added_keys: - if 'total_mass' not in converted_parameters.keys() and 'total_mass' not in added_keys: - converted_parameters['total_mass'] = chirp_mass_and_mass_ratio_to_total_mass( - converted_parameters['chirp_mass'], converted_parameters['mass_ratio']) - if remove: - added_keys.append('chirp_mass') - converted_parameters['mass_1'], converted_parameters['mass_2'] = \ - total_mass_and_mass_ratio_to_component_masses(converted_parameters['mass_ratio'], - converted_parameters['total_mass']) - if remove: - added_keys.append('total_mass') - added_keys.append('mass_ratio') - added_keys.append('mass_1') - added_keys.append('mass_2') - elif 'total_mass' in converted_parameters.keys() and 'mass_ratio' not in added_keys: - if 'symmetric_mass_ratio' in converted_parameters.keys() and 'symmetric_mass_ratio' not in added_keys: - converted_parameters['mass_ratio'] = \ - symmetric_mass_ratio_to_mass_ratio(converted_parameters['symmetric_mass_ratio']) - if remove: - added_keys.append('symmetric_mass_ratio') - if 'mass_ratio' in converted_parameters.keys() and 'mass_ratio' not in added_keys: - converted_parameters['mass_1'], converted_parameters['mass_2'] = \ - total_mass_and_mass_ratio_to_component_masses( - converted_parameters['mass_ratio'], converted_parameters['total_mass']) - if remove: - added_keys.append('total_mass') - added_keys.append('mass_ratio') - added_keys.append('mass_1') - added_keys.append('mass_2') - - elif 'mass_1' in search_keys and 'mass_2' not in search_keys: - if 'chirp_mass' in converted_parameters.keys() and 'chirp_mass' not in added_keys: + original_keys = list(converted_parameters.keys()) + + if 'chirp_mass' in converted_parameters.keys(): + if 'total_mass' in converted_parameters.keys(): + converted_parameters['symmetric_mass_ratio'] =\ + chirp_mass_and_total_mass_to_symmetric_mass_ratio( + converted_parameters['chirp_mass'], + converted_parameters['total_mass']) + if 'symmetric_mass_ratio' in converted_parameters.keys(): + converted_parameters['mass_ratio'] =\ + symmetric_mass_ratio_to_mass_ratio( + converted_parameters['symmetric_mass_ratio']) + if 'total_mass' not in converted_parameters.keys(): + converted_parameters['total_mass'] =\ + chirp_mass_and_mass_ratio_to_total_mass( + converted_parameters['chirp_mass'], + converted_parameters['mass_ratio']) + converted_parameters['mass_1'], converted_parameters['mass_2'] = \ + total_mass_and_mass_ratio_to_component_masses( + converted_parameters['mass_ratio'], + converted_parameters['total_mass']) + elif 'total_mass' in converted_parameters.keys(): + if 'symmetric_mass_ratio' in converted_parameters.keys(): converted_parameters['mass_ratio'] = \ - mass_1_and_chirp_mass_to_mass_ratio(parameters['mass_1'], parameters['chirp_mass']) - temp = (parameters['chirp_mass'] / parameters['mass_1']) ** 5 - parameters['mass_ratio'] = ( - (2 * temp / 3 / ( - (51 * temp ** 2 - 12 * temp ** 3) ** 0.5 + 9 * temp)) ** (1 / 3) + - (((51 * temp ** 2 - 12 * temp ** 3) ** 0.5 + 9 * temp) / 9 / 2 ** 0.5) ** (1 / 3)) - if remove: - added_keys.append('chirp_mass') - elif 'symmetric_mass_ratio' in converted_parameters.keys() and 'symmetric_mass_ratio' not in added_keys: - converted_parameters['mass_ratio'] = symmetric_mass_ratio_to_mass_ratio(parameters['symmetric_mass_ratio']) - if remove: - added_keys.append('symmetric_mass_ratio') - if 'mass_ratio' in converted_parameters.keys() and 'mass_ratio' not in added_keys: - converted_parameters['mass_2'] = converted_parameters['mass_1'] * converted_parameters['mass_ratio'] - if remove: - added_keys.append('mass_ratio') - added_keys.append('mass_2') - elif 'total_mass' in converted_parameters.keys() and 'total_mass' not in added_keys: - converted_parameters['mass_2'] = parameters['total_mass'] - parameters['mass_1'] - if remove: - added_keys.append('total_mass') - added_keys.append('mass_2') + symmetric_mass_ratio_to_mass_ratio( + converted_parameters['symmetric_mass_ratio']) + if 'mass_ratio' in converted_parameters.keys(): + converted_parameters['mass_1'], converted_parameters['mass_2'] =\ + total_mass_and_mass_ratio_to_component_masses( + converted_parameters['mass_ratio'], + converted_parameters['total_mass']) + elif 'mass_1' in converted_parameters.keys(): + converted_parameters['mass_2'] =\ + converted_parameters['total_mass'] -\ + converted_parameters['mass_1'] + elif 'mass_2' in converted_parameters.keys(): + converted_parameters['mass_1'] = \ + converted_parameters['total_mass'] - \ + converted_parameters['mass_2'] + elif 'symmetric_mass_ratio' in converted_parameters.keys(): + converted_parameters['mass_ratio'] =\ + symmetric_mass_ratio_to_mass_ratio( + converted_parameters['symmetric_mass_ratio']) + if 'mass_1' in converted_parameters.keys(): + converted_parameters['mass_2'] =\ + converted_parameters['mass_1'] *\ + converted_parameters['mass_ratio'] + elif 'mass_2' in converted_parameters.keys(): + converted_parameters['mass_1'] =\ + converted_parameters['mass_2'] /\ + converted_parameters['mass_ratio'] + elif 'mass_ratio' in converted_parameters.keys(): + if 'mass_1' in converted_parameters.keys(): + converted_parameters['mass_2'] =\ + converted_parameters['mass_1'] *\ + converted_parameters['mass_ratio'] + if 'mass_2' in converted_parameters.keys(): + converted_parameters['mass_1'] = \ + converted_parameters['mass_2'] /\ + converted_parameters['mass_ratio'] for angle in ['tilt_1', 'tilt_2', 'iota']: cos_angle = str('cos_' + angle) if cos_angle in converted_parameters.keys(): - converted_parameters[angle] = np.arccos(converted_parameters[cos_angle]) - added_keys.append(angle) - - if 'luminosity_distance' not in search_keys: - if 'redshift' in converted_parameters.keys(): - converted_parameters['luminosity_distance'] = redshift_to_luminosity_distance(parameters['redshift']) - if remove: - added_keys.append('redshift') - elif 'comoving_distance' in converted_parameters.keys(): - converted_parameters['luminosity_distance'] = \ - comoving_distance_to_luminosity_distance(parameters['comoving_distance']) - if remove: - added_keys.append('comoving_distance') - - added_keys = [key for key in added_keys if key not in search_keys] + converted_parameters[angle] =\ + np.arccos(converted_parameters[cos_angle]) + + if 'redshift' in converted_parameters.keys(): + converted_parameters['luminosity_distance'] =\ + redshift_to_luminosity_distance(parameters['redshift']) + elif 'comoving_distance' in converted_parameters.keys(): + converted_parameters['luminosity_distance'] = \ + comoving_distance_to_luminosity_distance( + parameters['comoving_distance']) + + added_keys = [key for key in converted_parameters.keys() + if key not in original_keys] return converted_parameters, added_keys -def convert_to_lal_binary_neutron_star_parameters(parameters, search_keys, remove=True): +def convert_to_lal_binary_neutron_star_parameters(parameters): """ Convert parameters we have into parameters we need. @@ -189,10 +174,6 @@ def convert_to_lal_binary_neutron_star_parameters(parameters, search_keys, remov ---------- parameters: dict dictionary of parameter values to convert into the required parameters - search_keys: list - parameters which are needed for the waveform generation - remove: bool, optional - Whether or not to remove the extra key, necessary for sampling, default=True. Return ------ @@ -202,36 +183,34 @@ def convert_to_lal_binary_neutron_star_parameters(parameters, search_keys, remov keys which are added to parameters during function call """ converted_parameters = parameters.copy() - converted_parameters, added_keys = convert_to_lal_binary_black_hole_parameters( - converted_parameters, search_keys, remove=remove) - - if 'lambda_1' not in search_keys and 'lambda_2' not in search_keys: - if 'delta_lambda' in converted_parameters.keys(): - converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\ - lambda_tilde_delta_lambda_to_lambda_1_lambda_2( - converted_parameters['lambda_tilde'], parameters['delta_lambda'], - converted_parameters['mass_1'], converted_parameters['mass_2']) - added_keys.append('lambda_1') - added_keys.append('lambda_2') - elif 'lambda_tilde' in converted_parameters.keys(): - converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\ - lambda_tilde_to_lambda_1_lambda_2( - converted_parameters['lambda_tilde'], - converted_parameters['mass_1'], converted_parameters['mass_2']) - added_keys.append('lambda_1') - added_keys.append('lambda_2') + original_keys = list(converted_parameters.keys()) + converted_parameters, added_keys =\ + convert_to_lal_binary_black_hole_parameters(converted_parameters) + + if 'delta_lambda' in converted_parameters.keys(): + converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\ + lambda_tilde_delta_lambda_to_lambda_1_lambda_2( + converted_parameters['lambda_tilde'], + parameters['delta_lambda'], converted_parameters['mass_1'], + converted_parameters['mass_2']) + elif 'lambda_tilde' in converted_parameters.keys(): + converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\ + lambda_tilde_to_lambda_1_lambda_2( + converted_parameters['lambda_tilde'], + converted_parameters['mass_1'], converted_parameters['mass_2']) if 'lambda_2' not in converted_parameters.keys(): converted_parameters['lambda_2'] =\ converted_parameters['lambda_1']\ * converted_parameters['mass_1']**5\ / converted_parameters['mass_2']**5 - added_keys.append('lambda_2') elif converted_parameters['lambda_2'] is None: converted_parameters['lambda_2'] =\ converted_parameters['lambda_1']\ * converted_parameters['mass_1']**5\ / converted_parameters['mass_2']**5 - added_keys.append('lambda_2') + + added_keys = [key for key in converted_parameters.keys() + if key not in original_keys] return converted_parameters, added_keys @@ -338,7 +317,7 @@ def component_masses_to_chirp_mass(mass_1, mass_2): Chirp mass of the binary """ - return (mass_1 * mass_2) ** 0.6 / (component_masses_to_total_mass(mass_1, mass_2)) ** 0.2 + return (mass_1 * mass_2) ** 0.6 / (mass_1 + mass_2) ** 0.2 def component_masses_to_total_mass(mass_1, mass_2): @@ -420,8 +399,10 @@ def mass_1_and_chirp_mass_to_mass_ratio(mass_1, chirp_mass): Mass ratio of the binary """ temp = (chirp_mass / mass_1) ** 5 - mass_ratio = (2 / 3 / (3 ** 0.5 * (27 * temp ** 2 - 4 * temp ** 3) ** 0.5 + 9 * temp)) ** (1 / 3) * temp + \ - ((3 ** 0.5 * (27 * temp ** 2 - 4 * temp ** 3) ** 0.5 + 9 * temp) / (2 * 3 ** 2)) ** (1 / 3) + mass_ratio = (2 / 3 / (3 ** 0.5 * (27 * temp ** 2 - 4 * temp ** 3) ** 0.5 + + 9 * temp)) ** (1 / 3) * temp + \ + ((3 ** 0.5 * (27 * temp ** 2 - 4 * temp ** 3) ** 0.5 + + 9 * temp) / (2 * 3 ** 2)) ** (1 / 3) return mass_ratio @@ -505,25 +486,32 @@ def lambda_tilde_to_lambda_1_lambda_2( def generate_all_bbh_parameters(sample, likelihood=None, priors=None): """ - From either a single sample or a set of samples fill in all missing BBH parameters, in place. + From either a single sample or a set of samples fill in all missing + BBH parameters, in place. Parameters ---------- sample: dict or pandas.DataFrame - Samples to fill in with extra parameters, this may be either an injection or posterior samples. - likelihood: tupak.gw.likelihood.GravitationalWaveTr, optional - GravitationalWaveTransient used for sampling, used for waveform and likelihood.interferometers. + Samples to fill in with extra parameters, this may be either an + injection or posterior samples. + likelihood: tupak.gw.likelihood.GravitationalWaveTransient, optional + GravitationalWaveTransient used for sampling, used for waveform and + likelihood.interferometers. priors: dict, optional Dictionary of prior objects, used to fill in non-sampled parameters. """ output_sample = sample.copy() if likelihood is not None: - output_sample['reference_frequency'] = likelihood.waveform_generator.waveform_arguments['reference_frequency'] - output_sample['waveform_approximant'] = likelihood.waveform_generator.waveform_arguments['waveform_approximant'] + output_sample['reference_frequency'] =\ + likelihood.waveform_generator.waveform_arguments[ + 'reference_frequency'] + output_sample['waveform_approximant'] =\ + likelihood.waveform_generator.waveform_arguments[ + 'waveform_approximant'] output_sample = fill_from_fixed_priors(output_sample, priors) - output_sample, _ = convert_to_lal_binary_black_hole_parameters( - output_sample, [key for key in output_sample.keys()], remove=False) + output_sample, _ =\ + convert_to_lal_binary_black_hole_parameters(output_sample) output_sample = generate_non_standard_parameters(output_sample) output_sample = generate_component_spins(output_sample) compute_snrs(output_sample, likelihood) @@ -557,7 +545,8 @@ def generate_non_standard_parameters(sample): Add the known non-standard parameters to the data frame/dictionary. We add: - chirp mass, total mass, symmetric mass ratio, mass ratio, cos tilt 1, cos tilt 2, cos iota + chirp mass, total mass, symmetric mass ratio, mass ratio, cos tilt 1, + cos tilt 2, cos iota Parameters ---------- @@ -570,16 +559,22 @@ def generate_non_standard_parameters(sample): """ output_sample = sample.copy() - output_sample['chirp_mass'] = component_masses_to_chirp_mass(sample['mass_1'], sample['mass_2']) - output_sample['total_mass'] = component_masses_to_total_mass(sample['mass_1'], sample['mass_2']) - output_sample['symmetric_mass_ratio'] = component_masses_to_symmetric_mass_ratio(sample['mass_1'], sample['mass_2']) - output_sample['mass_ratio'] = component_masses_to_mass_ratio(sample['mass_1'], sample['mass_2']) + output_sample['chirp_mass'] =\ + component_masses_to_chirp_mass(sample['mass_1'], sample['mass_2']) + output_sample['total_mass'] =\ + component_masses_to_total_mass(sample['mass_1'], sample['mass_2']) + output_sample['symmetric_mass_ratio'] =\ + component_masses_to_symmetric_mass_ratio(sample['mass_1'], + sample['mass_2']) + output_sample['mass_ratio'] =\ + component_masses_to_mass_ratio(sample['mass_1'], sample['mass_2']) output_sample['cos_tilt_1'] = np.cos(output_sample['tilt_1']) output_sample['cos_tilt_2'] = np.cos(output_sample['tilt_2']) output_sample['cos_iota'] = np.cos(output_sample['iota']) - output_sample['redshift'] = luminosity_distance_to_redshift(sample['luminosity_distance']) + output_sample['redshift'] =\ + luminosity_distance_to_redshift(sample['luminosity_distance']) return output_sample @@ -592,7 +587,8 @@ def generate_component_spins(sample): Parameters ---------- sample: A dictionary with the necessary spin conversion parameters: - 'iota', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1', 'mass_2', 'reference_frequency', 'phase' + 'iota', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1', + 'mass_2', 'reference_frequency', 'phase' Returns ------- @@ -600,45 +596,59 @@ def generate_component_spins(sample): """ output_sample = sample.copy() - spin_conversion_parameters = ['iota', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1', - 'mass_2', 'reference_frequency', 'phase'] + spin_conversion_parameters =\ + ['iota', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1', + 'mass_2', 'reference_frequency', 'phase'] if all(key in output_sample.keys() for key in spin_conversion_parameters)\ and isinstance(output_sample, dict): - output_sample['iota'], output_sample['spin_1x'], output_sample['spin_1y'], output_sample['spin_1z'], \ - output_sample['spin_2x'], output_sample['spin_2y'], output_sample['spin_2z'] = \ + output_sample['iota'], output_sample['spin_1x'],\ + output_sample['spin_1y'], output_sample['spin_1z'], \ + output_sample['spin_2x'], output_sample['spin_2y'],\ + output_sample['spin_2z'] =\ lalsim.SimInspiralTransformPrecessingNewInitialConditions( output_sample['iota'], output_sample['phi_jl'], output_sample['tilt_1'], output_sample['tilt_2'], - output_sample['phi_12'], output_sample['a_1'], output_sample['a_2'], + output_sample['phi_12'], output_sample['a_1'], + output_sample['a_2'], output_sample['mass_1'] * tupak.core.utils.solar_mass, output_sample['mass_2'] * tupak.core.utils.solar_mass, output_sample['reference_frequency'], output_sample['phase']) - output_sample['phi_1'] = np.arctan(output_sample['spin_1y'] / output_sample['spin_1x']) - output_sample['phi_2'] = np.arctan(output_sample['spin_2y'] / output_sample['spin_2x']) + output_sample['phi_1'] =\ + np.arctan(output_sample['spin_1y'] / output_sample['spin_1x']) + output_sample['phi_2'] =\ + np.arctan(output_sample['spin_2y'] / output_sample['spin_2x']) elif all(key in output_sample.keys() for key in spin_conversion_parameters)\ and isinstance(output_sample, pd.DataFrame): logger.debug('Extracting component spins.') - new_spin_parameters = ['spin_1x', 'spin_1y', 'spin_1z', 'spin_2x', 'spin_2y', 'spin_2z'] - new_spins = {name: np.zeros(len(output_sample)) for name in new_spin_parameters} + new_spin_parameters =\ + ['spin_1x', 'spin_1y', 'spin_1z', 'spin_2x', 'spin_2y', 'spin_2z'] + new_spins =\ + {name: np.zeros(len(output_sample)) for name in new_spin_parameters} for ii in range(len(output_sample)): - new_spins['iota'], new_spins['spin_1x'][ii], new_spins['spin_1y'][ii], new_spins['spin_1z'][ii], \ - new_spins['spin_2x'][ii], new_spins['spin_2y'][ii], new_spins['spin_2z'][ii] = \ + new_spins['iota'], new_spins['spin_1x'][ii],\ + new_spins['spin_1y'][ii], new_spins['spin_1z'][ii], \ + new_spins['spin_2x'][ii], new_spins['spin_2y'][ii],\ + new_spins['spin_2z'][ii] =\ lalsim.SimInspiralTransformPrecessingNewInitialConditions( output_sample['iota'][ii], output_sample['phi_jl'][ii], output_sample['tilt_1'][ii], output_sample['tilt_2'][ii], - output_sample['phi_12'][ii], output_sample['a_1'][ii], output_sample['a_2'][ii], + output_sample['phi_12'][ii], output_sample['a_1'][ii], + output_sample['a_2'][ii], output_sample['mass_1'][ii] * tupak.core.utils.solar_mass, output_sample['mass_2'][ii] * tupak.core.utils.solar_mass, - output_sample['reference_frequency'][ii], output_sample['phase'][ii]) + output_sample['reference_frequency'][ii], + output_sample['phase'][ii]) for name in new_spin_parameters: output_sample[name] = new_spins[name] - output_sample['phi_1'] = np.arctan(output_sample['spin_1y'] / output_sample['spin_1x']) - output_sample['phi_2'] = np.arctan(output_sample['spin_2y'] / output_sample['spin_2x']) + output_sample['phi_1'] =\ + np.arctan(output_sample['spin_1y'] / output_sample['spin_1x']) + output_sample['phi_2'] =\ + np.arctan(output_sample['spin_2y'] / output_sample['spin_2x']) else: logger.warning("Component spin extraction failed.") @@ -647,7 +657,9 @@ def generate_component_spins(sample): def compute_snrs(sample, likelihood): - """Compute the optimal and matched filter snrs of all posterior samples and print it out. + """ + Compute the optimal and matched filter snrs of all posterior samples + and print it out. Parameters ---------- @@ -660,41 +672,51 @@ def compute_snrs(sample, likelihood): temp_sample = sample if likelihood is not None: if isinstance(temp_sample, dict): + temp = dict() for key in likelihood.waveform_generator.parameters.keys(): - likelihood.waveform_generator.parameters[key] = temp_sample[key] - signal_polarizations = likelihood.waveform_generator.frequency_domain_strain() - for interferometer in likelihood.interferometers: - signal = interferometer.get_detector_response(signal_polarizations, - likelihood.waveform_generator.parameters) - sample['{}_matched_filter_snr'.format(interferometer.name)] = \ - interferometer.matched_filter_snr_squared(signal=signal) ** 0.5 - sample['{}_optimal_snr'.format(interferometer.name)] = \ - interferometer.optimal_snr_squared(signal=signal) ** 0.5 + # likelihood.waveform_generator.parameters[key] = temp_sample[key] + temp[key] = temp_sample[key] + signal_polarizations =\ + likelihood.waveform_generator.frequency_domain_strain(temp) + for ifo in likelihood.interferometers: + signal = ifo.get_detector_response( + signal_polarizations, + likelihood.waveform_generator.parameters) + sample['{}_matched_filter_snr'.format(ifo.name)] =\ + ifo.matched_filter_snr_squared(signal=signal) ** 0.5 + sample['{}_optimal_snr'.format(ifo.name)] = \ + ifo.optimal_snr_squared(signal=signal) ** 0.5 else: - logger.info('Computing SNRs for every sample, this may take some time.') + logger.info( + 'Computing SNRs for every sample, this may take some time.') all_interferometers = likelihood.interferometers - matched_filter_snrs = {interferometer.name: [] for interferometer in all_interferometers} - optimal_snrs = {interferometer.name: [] for interferometer in all_interferometers} + matched_filter_snrs = {ifo.name: [] for ifo in all_interferometers} + optimal_snrs = {ifo.name: [] for ifo in all_interferometers} for ii in range(len(temp_sample)): - for key in set(temp_sample.keys()).intersection(likelihood.waveform_generator.parameters.keys()): - likelihood.waveform_generator.parameters[key] = temp_sample[key][ii] - if likelihood.waveform_generator.non_standard_sampling_parameter_keys is not None: - for key in likelihood.waveform_generator.non_standard_sampling_parameter_keys: - likelihood.waveform_generator.parameters[key] = temp_sample[key][ii] - signal_polarizations = likelihood.waveform_generator.frequency_domain_strain() - for interferometer in all_interferometers: - signal = interferometer.get_detector_response(signal_polarizations, - likelihood.waveform_generator.parameters) - matched_filter_snrs[interferometer.name]\ - .append(interferometer.matched_filter_snr_squared(signal=signal) ** 0.5) - optimal_snrs[interferometer.name].append(interferometer.optimal_snr_squared(signal=signal) ** 0.5) - - for interferometer in likelihood.interferometers: - sample['{}_matched_filter_snr'.format(interferometer.name)] = matched_filter_snrs[interferometer.name] - sample['{}_optimal_snr'.format(interferometer.name)] = optimal_snrs[interferometer.name] + temp = dict() + for key in set(temp_sample.keys()).intersection( + likelihood.waveform_generator.parameters.keys()): + # likelihood.waveform_generator.parameters[key] =\ + # temp_sample[key][ii] + temp[key] = temp_sample[key][ii] + signal_polarizations =\ + likelihood.waveform_generator.frequency_domain_strain(temp) + for ifo in all_interferometers: + signal = ifo.get_detector_response( + signal_polarizations, + likelihood.waveform_generator.parameters) + matched_filter_snrs[ifo.name].append( + ifo.matched_filter_snr_squared(signal=signal) ** 0.5) + optimal_snrs[ifo.name].append( + ifo.optimal_snr_squared(signal=signal) ** 0.5) + + for ifo in likelihood.interferometers: + sample['{}_matched_filter_snr'.format(ifo.name)] =\ + matched_filter_snrs[ifo.name] + sample['{}_optimal_snr'.format(ifo.name)] =\ + optimal_snrs[ifo.name] likelihood.interferometers = all_interferometers - print([interferometer.name for interferometer in likelihood.interferometers]) else: logger.debug('Not computing SNRs.') diff --git a/tupak/gw/detector.py b/tupak/gw/detector.py index b571d7cc55e8912e3948db7d8dba62ad338e9798..924e8a8782149439c956381d431edd0114ac80dd 100644 --- a/tupak/gw/detector.py +++ b/tupak/gw/detector.py @@ -107,8 +107,8 @@ class InterferometerList(list): """ if injection_polarizations is None: if waveform_generator is not None: - waveform_generator.parameters = parameters - injection_polarizations = waveform_generator.frequency_domain_strain() + injection_polarizations =\ + waveform_generator.frequency_domain_strain(parameters) else: raise ValueError( "inject_signal needs one of waveform_generator or " @@ -1241,8 +1241,8 @@ class Interferometer(object): if injection_polarizations is None: if waveform_generator is not None: - waveform_generator.parameters = parameters - injection_polarizations = waveform_generator.frequency_domain_strain() + injection_polarizations =\ + waveform_generator.frequency_domain_strain(parameters) else: raise ValueError( "inject_signal needs one of waveform_generator or " diff --git a/tupak/gw/likelihood.py b/tupak/gw/likelihood.py index ec89d11f08655f85b845b8c24135b1fb47a1a4d4..3d7a77197c89cde2caa64042dd663ab7567ddb1d 100644 --- a/tupak/gw/likelihood.py +++ b/tupak/gw/likelihood.py @@ -56,7 +56,7 @@ class GravitationalWaveTransient(likelihood.Likelihood): phase_marginalization=False, prior=None): self.waveform_generator = waveform_generator - likelihood.Likelihood.__init__(self, waveform_generator.parameters) + likelihood.Likelihood.__init__(self, dict()) self.interferometers = tupak.gw.detector.InterferometerList(interferometers) self.time_marginalization = time_marginalization self.distance_marginalization = distance_marginalization @@ -135,18 +135,6 @@ class GravitationalWaveTransient(likelihood.Likelihood): else: self.__prior = None - @property - def non_standard_sampling_parameter_keys(self): - return self.waveform_generator.non_standard_sampling_parameter_keys - - @property - def parameters(self): - return self.waveform_generator.parameters - - @parameters.setter - def parameters(self, parameters): - self.waveform_generator.parameters = parameters - def noise_log_likelihood(self): log_l = 0 for interferometer in self.interferometers: @@ -159,7 +147,7 @@ class GravitationalWaveTransient(likelihood.Likelihood): def log_likelihood_ratio(self): waveform_polarizations =\ - self.waveform_generator.frequency_domain_strain() + self.waveform_generator.frequency_domain_strain(self.parameters) if waveform_polarizations is None: return np.nan_to_num(-np.inf) @@ -171,7 +159,7 @@ class GravitationalWaveTransient(likelihood.Likelihood): dtype=np.complex128) for interferometer in self.interferometers: signal_ifo = interferometer.get_detector_response( - waveform_polarizations, self.waveform_generator.parameters) + waveform_polarizations, self.parameters) matched_filter_snr_squared += interferometer.matched_filter_snr_squared(signal=signal_ifo) optimal_snr_squared += interferometer.optimal_snr_squared(signal=signal_ifo) @@ -220,10 +208,10 @@ class GravitationalWaveTransient(likelihood.Likelihood): def _setup_rho(self, matched_filter_snr_squared, optimal_snr_squared): rho_opt_ref = (optimal_snr_squared.real * - self.waveform_generator.parameters['luminosity_distance'] ** 2 / + self.parameters['luminosity_distance'] ** 2 / self._ref_dist ** 2.) rho_mf_ref = matched_filter_snr_squared * \ - self.waveform_generator.parameters['luminosity_distance'] / self._ref_dist + self.parameters['luminosity_distance'] / self._ref_dist return rho_mf_ref, rho_opt_ref def log_likelihood(self): @@ -310,7 +298,7 @@ class BasicGravitationalWaveTransient(likelihood.Likelihood): given some set of parameters """ - likelihood.Likelihood.__init__(self, waveform_generator.parameters) + likelihood.Likelihood.__init__(self, dict()) self.interferometers = interferometers self.waveform_generator = waveform_generator @@ -342,7 +330,9 @@ class BasicGravitationalWaveTransient(likelihood.Likelihood): """ log_l = 0 - waveform_polarizations = self.waveform_generator.frequency_domain_strain() + waveform_polarizations =\ + self.waveform_generator.frequency_domain_strain( + self.parameters.copy()) if waveform_polarizations is None: return np.nan_to_num(-np.inf) for interferometer in self.interferometers: @@ -367,7 +357,7 @@ class BasicGravitationalWaveTransient(likelihood.Likelihood): """ signal_ifo = interferometer.get_detector_response( - waveform_polarizations, self.waveform_generator.parameters) + waveform_polarizations, self.parameters) log_l = - 2. / self.waveform_generator.duration * np.vdot( interferometer.frequency_domain_strain - signal_ifo, @@ -394,5 +384,5 @@ def get_binary_black_hole_likelihood(interferometers): waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( duration=interferometers.duration, sampling_frequency=interferometers.sampling_frequency, frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, - parameters={'waveform_approximant': 'IMRPhenomPv2', 'reference_frequency': 50}) + waveform_arguments={'waveform_approximant': 'IMRPhenomPv2', 'reference_frequency': 50}) return tupak.gw.likelihood.GravitationalWaveTransient(interferometers, waveform_generator) diff --git a/tupak/gw/waveform_generator.py b/tupak/gw/waveform_generator.py index b5e08cc15548c53ba661ab3393db8925509f281c..0134ffbb5f0c7ae6caeffbb830a4308ff5086c50 100644 --- a/tupak/gw/waveform_generator.py +++ b/tupak/gw/waveform_generator.py @@ -7,7 +7,6 @@ class WaveformGenerator(object): def __init__(self, duration=None, sampling_frequency=None, start_time=0, frequency_domain_source_model=None, time_domain_source_model=None, parameters=None, parameter_conversion=None, - non_standard_sampling_parameter_keys=None, waveform_arguments=None): """ A waveform generator @@ -33,8 +32,6 @@ class WaveformGenerator(object): Function to convert from sampled parameters to parameters of the waveform generator. Default value is the identity, i.e. it leaves the parameters unaffected. - non_standard_sampling_parameter_keys: list, optional - List of parameter name for *non-standard* sampling parameters. waveform_arguments: dict, optional A dictionary of fixed keyword arguments to pass to either `frequency_domain_source_model` or `time_domain_source_model`. @@ -49,25 +46,21 @@ class WaveformGenerator(object): self.start_time = start_time self.frequency_domain_source_model = frequency_domain_source_model self.time_domain_source_model = time_domain_source_model - self.__parameters_from_source_model() + self.source_parameter_keys = self.__parameters_from_source_model() self.duration = duration self.sampling_frequency = sampling_frequency if parameter_conversion is None: - self.parameter_conversion = lambda params, search_keys: (params, []) + self.parameter_conversion = lambda params: (params, []) else: self.parameter_conversion = parameter_conversion - self.non_standard_sampling_parameter_keys = non_standard_sampling_parameter_keys - self.parameters = parameters if waveform_arguments is not None: self.waveform_arguments = waveform_arguments else: self.waveform_arguments = dict() + if isinstance(parameters, dict): + self.parameters = parameters self.__frequency_array_updated = False self.__time_array_updated = False - self.__full_source_model_keyword_arguments = {} - self.__full_source_model_keyword_arguments.update(self.waveform_arguments) - self.__full_source_model_keyword_arguments.update(self.parameters) - self.__added_keys = [] def __repr__(self): if self.frequency_domain_source_model is not None: @@ -85,17 +78,24 @@ class WaveformGenerator(object): return self.__class__.__name__ + '(duration={}, sampling_frequency={}, start_time={}, ' \ 'frequency_domain_source_model={}, time_domain_source_model={}, ' \ - 'parameters={}, parameter_conversion={}, ' \ - 'non_standard_sampling_parameter_keys={}, waveform_arguments={})'\ - .format(self.duration, self.sampling_frequency, self.start_time, fdsm_name, tdsm_name, self.parameters, - param_conv_name, self.non_standard_sampling_parameter_keys, self.waveform_arguments) + 'parameter_conversion={}, ' \ + 'waveform_arguments={})'\ + .format(self.duration, self.sampling_frequency, self.start_time, fdsm_name, tdsm_name, + param_conv_name, self.waveform_arguments) - def frequency_domain_strain(self): - """ Rapper to source_model. + def frequency_domain_strain(self, parameters=None): + """ Wrapper to source_model. Converts self.parameters with self.parameter_conversion before handing it off to the source model. Automatically refers to the time_domain_source model via NFFT if no frequency_domain_source_model is given. + Parameters + ---------- + parameters: dict, optional + Parameters to evaluate the waveform for, this overwrites + `self.parameters`. + If not provided will fall back to `self.parameters`. + Returns ------- array_like: The frequency domain strain for the given set of parameters @@ -107,17 +107,25 @@ class WaveformGenerator(object): """ return self._calculate_strain(model=self.frequency_domain_source_model, model_data_points=self.frequency_array, + parameters=parameters, transformation_function=utils.nfft, transformed_model=self.time_domain_source_model, transformed_model_data_points=self.time_array) - def time_domain_strain(self): - """ Rapper to source_model. + def time_domain_strain(self, parameters=None): + """ Wrapper to source_model. Converts self.parameters with self.parameter_conversion before handing it off to the source model. Automatically refers to the frequency_domain_source model via INFFT if no frequency_domain_source_model is given. + Parameters + ---------- + parameters: dict, optional + Parameters to evaluate the waveform for, this overwrites + `self.parameters`. + If not provided will fall back to `self.parameters`. + Returns ------- array_like: The time domain strain for the given set of parameters @@ -129,13 +137,15 @@ class WaveformGenerator(object): """ return self._calculate_strain(model=self.time_domain_source_model, model_data_points=self.time_array, + parameters=parameters, transformation_function=utils.infft, transformed_model=self.frequency_domain_source_model, transformed_model_data_points=self.frequency_array) def _calculate_strain(self, model, model_data_points, transformation_function, transformed_model, - transformed_model_data_points): - self._apply_parameter_conversion() + transformed_model_data_points, parameters): + if parameters is not None: + self.parameters = parameters if model is not None: model_strain = self._strain_from_model(model_data_points, model) elif transformed_model is not None: @@ -143,16 +153,10 @@ class WaveformGenerator(object): transformation_function) else: raise RuntimeError("No source model given") - self._remove_added_keys() return model_strain - def _apply_parameter_conversion(self): - self.parameters, self.__added_keys = self.parameter_conversion(self.parameters, - self.non_standard_sampling_parameter_keys) - self.__full_source_model_keyword_arguments.update(self.parameters) - def _strain_from_model(self, model_data_points, model): - return model(model_data_points, **self.__full_source_model_keyword_arguments) + return model(model_data_points, **self.parameters) def _strain_from_transformed_model(self, transformed_model_data_points, transformed_model, transformation_function): transformed_model_strain = self._strain_from_model(transformed_model_data_points, transformed_model) @@ -169,10 +173,6 @@ class WaveformGenerator(object): model_strain[key] = transformation_function(transformed_model_strain[key], self.sampling_frequency) return model_strain - def _remove_added_keys(self): - for key in self.__added_keys: - self.parameters.pop(key) - @property def frequency_array(self): """ Frequency array for the waveforms. Automatically updates if sampling_frequency or duration are updated. @@ -228,15 +228,44 @@ class WaveformGenerator(object): @parameters.setter def parameters(self, parameters): - if isinstance(parameters, dict): - for key in parameters.keys(): - self.__parameters[key] = parameters[key] + """ + Set parameters, this applies the conversion function and then removes + any parameters which aren't required by the source function. + + (set.symmetric_difference is the opposite of set.intersection) + + Parameters + ---------- + parameters: dict + Input parameter dictionary, this is copied, passed to the conversion + function and has self.waveform_arguments added to it. + """ + if not isinstance(parameters, dict): + raise TypeError('"parameters" must be a dictionary.') + new_parameters = parameters.copy() + new_parameters, _ = self.parameter_conversion(new_parameters) + for key in self.source_parameter_keys.symmetric_difference( + new_parameters): + new_parameters.pop(key) + self.__parameters = new_parameters + self.__parameters.update(self.waveform_arguments) def __parameters_from_source_model(self): + """ + Infer the named arguments of the source model. + + Returns + ------- + set: The names of the arguments of the source model. + """ if self.frequency_domain_source_model is not None: - self.__parameters = dict.fromkeys(utils.infer_parameters_from_function(self.frequency_domain_source_model)) + model = self.frequency_domain_source_model elif self.time_domain_source_model is not None: - self.__parameters = dict.fromkeys(utils.infer_parameters_from_function(self.time_domain_source_model)) + model = self.time_domain_source_model + else: + raise AttributeError('Either time or frequency domain source ' + 'model must be provided.') + return set(utils.infer_parameters_from_function(model)) @property def duration(self):