diff --git a/examples/injection_examples/basic_tutorial.py b/examples/injection_examples/basic_tutorial.py index c3548b4d41ad42bd5584c2bf863f00ffd79b385e..26923d0149a03dc39662a9242bf1c1eb9f47f333 100644 --- a/examples/injection_examples/basic_tutorial.py +++ b/examples/injection_examples/basic_tutorial.py @@ -13,7 +13,7 @@ import tupak # Set the duration and sampling frequency of the data segment that we're going to inject the signal into -time_duration = 4. +duration = 4. sampling_frequency = 2048. # Specify the output directory and the name of the simulation. @@ -36,7 +36,7 @@ waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', reference_frequency=50.) # Create the waveform_generator using a LAL BinaryBlackHole source function -waveform_generator = tupak.WaveformGenerator(time_duration=time_duration, +waveform_generator = tupak.WaveformGenerator(duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, parameters=injection_parameters, @@ -46,7 +46,7 @@ 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, + name, injection_polarizations=hf_signal, injection_parameters=injection_parameters, duration=duration, sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1']] # Set up prior, which is a dictionary diff --git a/examples/injection_examples/change_sampled_parameters.py b/examples/injection_examples/change_sampled_parameters.py index a563a32882f7b69b6f924fecd7a8e9cbde92f4d3..9f95dce802453bbd838f760db3d164b22a409b17 100644 --- a/examples/injection_examples/change_sampled_parameters.py +++ b/examples/injection_examples/change_sampled_parameters.py @@ -12,7 +12,7 @@ import numpy as np tupak.core.utils.setup_logger(log_level="info") -time_duration = 4. +duration = 4. sampling_frequency = 2048. outdir = 'outdir' @@ -27,7 +27,7 @@ waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', # Create the waveform_generator using a LAL BinaryBlackHole source function waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( - sampling_frequency=sampling_frequency, time_duration=time_duration, + 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'], @@ -36,7 +36,7 @@ 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, time_duration=time_duration, + name, injection_polarizations=hf_signal, injection_parameters=injection_parameters, duration=duration, sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']] # Set up prior diff --git a/examples/injection_examples/create_your_own_source_model.py b/examples/injection_examples/create_your_own_source_model.py index a595695044d828ef7205ee13fc72eafdec043e29..604a513ef6bb9a70ba60a7db02e4d3dfecf6696f 100644 --- a/examples/injection_examples/create_your_own_source_model.py +++ b/examples/injection_examples/create_your_own_source_model.py @@ -10,7 +10,7 @@ import numpy as np outdir = 'outdir' label = 'create_your_own_source_model' sampling_frequency = 4096 -time_duration = 1 +duration = 1 # Here we define out source model - this is the sine-Gaussian model in the @@ -25,7 +25,7 @@ def sine_gaussian(f, A, f0, tau, phi0, geocent_time, ra, dec, psi): # We now define some parameters that we will inject and then a waveform generator injection_parameters = dict(A=1e-23, f0=100, tau=1, phi0=0, geocent_time=0, ra=0, dec=0, psi=0) -waveform_generator = tupak.gw.waveform_generator.WaveformGenerator(time_duration=time_duration, +waveform_generator = tupak.gw.waveform_generator.WaveformGenerator(duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=sine_gaussian, parameters=injection_parameters) @@ -34,7 +34,7 @@ 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, time_duration=time_duration, + injection_parameters=injection_parameters, duration=duration, sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']] 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 f89bb0f77c3ccbec9a9e39bd54bb23232ad7af06..be48be57ab9c23f8594da2751f03e5be586fb698 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 @@ -26,13 +26,13 @@ injection_parameters = dict(amplitude=5e-22, damping_time=0.1, frequency=50, phase=0, ra=0, dec=0, psi=0, geocent_time=0.) -time_duration = 0.5 +duration = 0.5 sampling_frequency = 2048 outdir='outdir' label='time_domain_source_model' # call the waveform_generator to create our waveform model. -waveform = tupak.gw.waveform_generator.WaveformGenerator(time_duration=time_duration, sampling_frequency=sampling_frequency, +waveform = tupak.gw.waveform_generator.WaveformGenerator(duration=duration, sampling_frequency=sampling_frequency, time_domain_source_model=time_domain_damped_sinusoid, parameters=injection_parameters) @@ -49,7 +49,7 @@ hf_signal = waveform.frequency_domain_strain() # 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, time_duration=time_duration, + injection_parameters=injection_parameters, duration=duration, sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1']] diff --git a/examples/injection_examples/how_to_specify_the_prior.py b/examples/injection_examples/how_to_specify_the_prior.py index d33fc753486bd01f8fbef8968a2cf407da54ecea..76b2884eacdf6ea6e1fd4f8f258bbb399045f90d 100644 --- a/examples/injection_examples/how_to_specify_the_prior.py +++ b/examples/injection_examples/how_to_specify_the_prior.py @@ -9,7 +9,7 @@ import numpy as np import tupak.gw.prior -time_duration = 4. +duration = 4. sampling_frequency = 2048. outdir = 'outdir' @@ -23,7 +23,7 @@ waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', reference_frequency=50.) # Create the waveform_generator using a LAL BinaryBlackHole source function -waveform_generator = tupak.WaveformGenerator(time_duration=time_duration, +waveform_generator = tupak.WaveformGenerator(duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, parameters=injection_parameters, @@ -32,7 +32,7 @@ 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, time_duration=time_duration, + name, injection_polarizations=hf_signal, injection_parameters=injection_parameters, duration=duration, sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']] # Set up prior diff --git a/examples/injection_examples/marginalized_likelihood.py b/examples/injection_examples/marginalized_likelihood.py index adf4cd9b4ffb43284e006e2d06d15a8c5e0d9caa..b2280bca42716cf03aca32913692e833c9dda666 100644 --- a/examples/injection_examples/marginalized_likelihood.py +++ b/examples/injection_examples/marginalized_likelihood.py @@ -8,7 +8,7 @@ import tupak import numpy as np -time_duration = 4. +duration = 4. sampling_frequency = 2048. outdir = 'outdir' @@ -23,14 +23,14 @@ waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', # Create the waveform_generator using a LAL BinaryBlackHole source function waveform_generator = tupak.WaveformGenerator( - time_duration=time_duration, sampling_frequency=sampling_frequency, + 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, time_duration=time_duration, + name, injection_polarizations=hf_signal, injection_parameters=injection_parameters, duration=duration, sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']] # Set up prior diff --git a/examples/injection_examples/sine_gaussian_example.py b/examples/injection_examples/sine_gaussian_example.py index adc1e9bbb1eb7dd17b147d32d04563a7c84f2c5d..8761cf6b388542a34717ec52c1c71c637d6fac1f 100644 --- a/examples/injection_examples/sine_gaussian_example.py +++ b/examples/injection_examples/sine_gaussian_example.py @@ -8,7 +8,7 @@ import tupak import numpy as np # Set the duration and sampling frequency of the data segment that we're going to inject the signal into -time_duration = 4. +duration = 4. sampling_frequency = 2048. # Specify the output directory and the name of the simulation. @@ -25,7 +25,7 @@ injection_parameters = dict(hrss = 1e-22, Q = 5.0, frequency = 200.0, ra = 1.375 geocent_time = 1126259642.413, psi= 2.659) # Create the waveform_generator using a sine Gaussian source function -waveform_generator = tupak.gw.waveform_generator.WaveformGenerator(time_duration=time_duration, +waveform_generator = tupak.gw.waveform_generator.WaveformGenerator(duration=duration, sampling_frequency=sampling_frequency, frequency_domain_source_model=tupak.gw.source.sinegaussian, parameters=injection_parameters) @@ -34,7 +34,7 @@ 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, + name, injection_polarizations=hf_signal, injection_parameters=injection_parameters, duration=duration, sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']] # Set up prior, which is a dictionary diff --git a/examples/open_data_examples/GW150914.py b/examples/open_data_examples/GW150914.py index f2f21b688c6691d09d50215e51e4fa4e4c9b6147..07a31db0209fdac49ecaa46ba7a2411b05d0c9bc 100644 --- a/examples/open_data_examples/GW150914.py +++ b/examples/open_data_examples/GW150914.py @@ -36,9 +36,7 @@ prior = tupak.gw.prior.BBHPriorSet(filename='GW150914.prior') # creates the frequency-domain strain. In this instance, we are using the # `lal_binary_black_hole model` source model. We also pass other parameters: # the waveform approximant and reference frequency. -waveform_generator = tupak.WaveformGenerator(time_duration=interferometers.duration, - sampling_frequency=interferometers.sampling_frequency, - frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, +waveform_generator = tupak.WaveformGenerator(frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, waveform_arguments={'waveform_approximant': 'IMRPhenomPv2', 'reference_frequency': 50}) diff --git a/test/waveform_generator_tests.py b/test/waveform_generator_tests.py index a13c284060d5cbf4bccd28765b08bbcced34336d..16c13c5dbb0113b38fa6f9df34f2ae05d10d84e5 100644 --- a/test/waveform_generator_tests.py +++ b/test/waveform_generator_tests.py @@ -40,8 +40,8 @@ class TestWaveformGeneratorInstantiationWithoutOptionalParameters(unittest.TestC del self.waveform_generator del self.simulation_parameters - def test_time_duration(self): - self.assertEqual(self.waveform_generator.time_duration, 1) + def test_duration(self): + self.assertEqual(self.waveform_generator.duration, 1) def test_sampling_frequency(self): self.assertEqual(self.waveform_generator.sampling_frequency, 4096) diff --git a/tupak/core/sampler.py b/tupak/core/sampler.py index 22e0797ec67c54623b119a2244a77ebe5f5eb3de..89d542dfff01a280ffc59b526754b717c09eb016 100644 --- a/tupak/core/sampler.py +++ b/tupak/core/sampler.py @@ -348,8 +348,15 @@ class Sampler(object): def _log_summary_for_sampler(self): """Print a summary of the sampler used and its kwargs""" if self.cached_result is None: + kwargs_print = self.kwargs.copy() + for k in kwargs_print: + if type(kwargs_print[k]) in (list, np.ndarray): + array_repr = np.array(kwargs_print[k]) + if array_repr.shape > 10: + kwargs_print[k] = ('array_like, shape={}' + .format(array_repr.shape)) logging.info("Using sampler {} with kwargs {}".format( - self.__class__.__name__, self.kwargs)) + self.__class__.__name__, kwargs_print)) class Nestle(Sampler): @@ -764,7 +771,17 @@ class Emcee(Sampler): sampler = emcee.EnsembleSampler( nwalkers=self.nwalkers, dim=self.ndim, lnpostfn=self.lnpostfn, a=a) - pos0 = [self.get_random_draw_from_prior() for i in range(self.nwalkers)] + + if 'pos0' in self.kwargs: + logging.debug("Using given initial positions for walkers") + pos0 = np.squeeze(self.kwargs['pos0']) + if pos0.shape != (self.ndim, self.nwalkers): + raise ValueError( + 'Input pos0 should be of shape ndim, nwalkers') + else: + logging.debug("Generating initial walker positions from prior") + pos0 = [self.get_random_draw_from_prior() + for i in range(self.nwalkers)] for result in tqdm( sampler.sample(pos0, iterations=self.nsteps), total=self.nsteps): diff --git a/tupak/gw/conversion.py b/tupak/gw/conversion.py index f08d3d0f59d185340bb0f2808ec90737a8a53a3d..39fcaec6cf5bc49044b031f6898c9b01307d32af 100644 --- a/tupak/gw/conversion.py +++ b/tupak/gw/conversion.py @@ -515,9 +515,9 @@ def compute_snrs(sample, likelihood): likelihood.waveform_generator.parameters) sample['{}_matched_filter_snr'.format(interferometer.name)] = \ tupak.gw.utils.matched_filter_snr_squared(signal, interferometer, - likelihood.waveform_generator.time_duration) ** 0.5 + likelihood.waveform_generator.duration) ** 0.5 sample['{}_optimal_snr'.format(interferometer.name)] = tupak.gw.utils.optimal_snr_squared( - signal, interferometer, likelihood.waveform_generator.time_duration) ** 0.5 + signal, interferometer, likelihood.waveform_generator.duration) ** 0.5 else: logging.info('Computing SNRs for every sample, this may take some time.') all_interferometers = likelihood.interferometers @@ -533,9 +533,9 @@ def compute_snrs(sample, likelihood): signal = interferometer.get_detector_response(signal_polarizations, likelihood.waveform_generator.parameters) matched_filter_snrs[interferometer.name].append(tupak.gw.utils.matched_filter_snr_squared( - signal, interferometer, likelihood.waveform_generator.time_duration) ** 0.5) + signal, interferometer, likelihood.waveform_generator.duration) ** 0.5) optimal_snrs[interferometer.name].append(tupak.gw.utils.optimal_snr_squared( - signal, interferometer, likelihood.waveform_generator.time_duration) ** 0.5) + signal, interferometer, likelihood.waveform_generator.duration) ** 0.5) for interferometer in likelihood.interferometers: sample['{}_matched_filter_snr'.format(interferometer.name)] = matched_filter_snrs[interferometer.name] diff --git a/tupak/gw/detector.py b/tupak/gw/detector.py index c6e78d3ca0ecdbf0667377443898f34b0e2eb676..ecbe133e2804f15f53f7446284e06075218cf643 100644 --- a/tupak/gw/detector.py +++ b/tupak/gw/detector.py @@ -539,11 +539,11 @@ class InterferometerStrainData(object): logging.debug( 'Setting data using noise realization from provided' 'power_spectal_density') - frequency_domain_strain, frequencies = \ + frequency_domain_strain, frequency_array = \ power_spectral_density.get_noise_realisation( self.sampling_frequency, self.duration) - if np.array_equal(frequencies, self.frequency_array): + if np.array_equal(frequency_array, self.frequency_array): self._frequency_domain_strain = frequency_domain_strain else: raise ValueError("Data frequencies do not match frequency_array") @@ -1107,10 +1107,10 @@ class Interferometer(object): start_time=self.strain_data.start_time) opt_snr = np.sqrt(tupak.gw.utils.optimal_snr_squared( signal=signal_ifo, interferometer=self, - time_duration=self.strain_data.duration).real) + duration=self.strain_data.duration).real) mf_snr = np.sqrt(tupak.gw.utils.matched_filter_snr_squared( signal=signal_ifo, interferometer=self, - time_duration=self.strain_data.duration).real) + duration=self.strain_data.duration).real) logging.info("Injected signal in {}:".format(self.name)) logging.info(" optimal SNR = {:.2f}".format(opt_snr)) @@ -1320,7 +1320,7 @@ class PowerSpectralDensity(object): Array representation of the ASD amplitude_spectral_density_file: str Name of the ASD file - frequencies: array_like + frequency_array: array_like Array containing the frequencies of the ASD/PSD values power_spectral_density: array_like Array representation of the PSD @@ -1333,7 +1333,7 @@ class PowerSpectralDensity(object): self.__power_spectral_density = None self.__amplitude_spectral_density = None - self.frequencies = [] + self.frequency_array = [] self.power_spectral_density_interpolated = None for key in kwargs: @@ -1417,16 +1417,16 @@ class PowerSpectralDensity(object): strain.low_pass_filter(filter_freq) f, psd = strain.create_power_spectral_density(fft_length=fft_length) - self.frequencies = f + self.frequency_array = f self.power_spectral_density = psd - def set_from_amplitude_spectral_density_array(self, frequencies, + def set_from_amplitude_spectral_density_array(self, frequency_array, asd_array): - self.frequencies = frequencies + self.frequency_array = frequency_array self.amplitude_spectral_density = asd_array - def set_from_power_spectral_density_array(self, frequencies, psd_array): - self.frequencies = frequencies + def set_from_power_spectral_density_array(self, frequency_array, psd_array): + self.frequency_array = frequency_array self.power_spectral_density = psd_array def set_from_aLIGO(self): @@ -1475,7 +1475,7 @@ class PowerSpectralDensity(object): os.path.dirname(__file__), 'noise_curves', self.amplitude_spectral_density_file) - self.frequencies, self.amplitude_spectral_density = np.genfromtxt( + self.frequency_array, self.amplitude_spectral_density = np.genfromtxt( self.amplitude_spectral_density_file).T def import_power_spectral_density(self): @@ -1490,13 +1490,13 @@ class PowerSpectralDensity(object): self.power_spectral_density_file = os.path.join( os.path.dirname(__file__), 'noise_curves', self.power_spectral_density_file) - self.frequencies, self.power_spectral_density = np.genfromtxt( + self.frequency_array, self.power_spectral_density = np.genfromtxt( self.power_spectral_density_file).T def _check_frequency_array_matches_density_array(self, density_array): """Check the provided frequency and spectral density arrays match.""" try: - self.frequencies - density_array + self.frequency_array - density_array except ValueError as e: raise(e, 'Provided spectral density does not match frequency array. Not updating.') @@ -1505,7 +1505,7 @@ class PowerSpectralDensity(object): for arbitrary frequency arrays. """ self.power_spectral_density_interpolated = interp1d( - self.frequencies, self.power_spectral_density, bounds_error=False, + self.frequency_array, self.power_spectral_density, bounds_error=False, fill_value=np.inf) def get_noise_realisation(self, sampling_frequency, duration): @@ -1528,7 +1528,7 @@ class PowerSpectralDensity(object): white_noise, frequencies = utils.create_white_noise(sampling_frequency, duration) interpolated_power_spectral_density = self.power_spectral_density_interpolated(frequencies) frequency_domain_strain = interpolated_power_spectral_density ** 0.5 * white_noise - out_of_bounds = (frequencies < min(self.frequencies)) | (frequencies > max(self.frequencies)) + out_of_bounds = (frequencies < min(self.frequency_array)) | (frequencies > max(self.frequency_array)) frequency_domain_strain[out_of_bounds] = 0 * (1 + 1j) return frequency_domain_strain, frequencies @@ -1588,7 +1588,7 @@ def load_interferometer(filename): def get_interferometer_with_open_data( - name, trigger_time, time_duration=4, start_time=None, roll_off=0.4, psd_offset=-1024, + name, trigger_time, duration=4, start_time=None, roll_off=0.4, psd_offset=-1024, psd_duration=100, cache=True, outdir='outdir', label=None, plot=True, filter_freq=None, raw_data_file=None, **kwargs): """ @@ -1602,7 +1602,7 @@ def get_interferometer_with_open_data( Detector name, e.g., 'H1'. trigger_time: float Trigger GPS time. - time_duration: float, optional + duration: float, optional The total time (in seconds) to analyse. Defaults to 4s. start_time: float, optional Beginning of the segment, if None, the trigger is placed 2s before the end @@ -1642,16 +1642,16 @@ def get_interferometer_with_open_data( utils.check_directory_exists_and_if_not_mkdir(outdir) if start_time is None: - start_time = trigger_time + 2 - time_duration + start_time = trigger_time + 2 - duration strain = InterferometerStrainData(roll_off=roll_off) strain.set_from_open_data( - name=name, start_time=start_time, duration=time_duration, + name=name, start_time=start_time, duration=duration, outdir=outdir, cache=cache, **kwargs) strain_psd = InterferometerStrainData(roll_off=roll_off) strain_psd.set_from_open_data( - name=name, start_time=start_time + time_duration + psd_offset, + name=name, start_time=start_time + duration + psd_offset, duration=psd_duration, outdir=outdir, cache=cache, **kwargs) # Low pass filter strain_psd.low_pass_filter(filter_freq) @@ -1661,7 +1661,7 @@ def get_interferometer_with_open_data( interferometer = get_empty_interferometer(name) interferometer.power_spectral_density = PowerSpectralDensity( - psd_array=psd_array, frequencies=psd_frequencies) + psd_array=psd_array, frequency_array=psd_frequencies) interferometer.strain_data = strain if plot: @@ -1672,7 +1672,7 @@ def get_interferometer_with_open_data( def get_interferometer_with_fake_noise_and_injection( name, injection_parameters, injection_polarizations=None, - waveform_generator=None, sampling_frequency=4096, time_duration=4, + waveform_generator=None, sampling_frequency=4096, duration=4, start_time=None, outdir='outdir', label=None, plot=True, save=True, zero_noise=False): """ @@ -1698,7 +1698,7 @@ def get_interferometer_with_fake_noise_and_injection( `injection_polarizations` is given, this will be ignored. sampling_frequency: float sampling frequency for data, should match injection signal - time_duration: float + duration: float length of data, should be the same as used for signal generation start_time: float Beginning of data segment, if None, injection is placed 2s before @@ -1723,17 +1723,17 @@ def get_interferometer_with_fake_noise_and_injection( utils.check_directory_exists_and_if_not_mkdir(outdir) if start_time is None: - start_time = injection_parameters['geocent_time'] + 2 - time_duration + start_time = injection_parameters['geocent_time'] + 2 - duration interferometer = get_empty_interferometer(name) interferometer.power_spectral_density.set_from_aLIGO() if zero_noise: interferometer.set_strain_data_from_zero_noise( - sampling_frequency=sampling_frequency, duration=time_duration, + sampling_frequency=sampling_frequency, duration=duration, start_time=start_time) else: interferometer.set_strain_data_from_power_spectral_density( - sampling_frequency=sampling_frequency, duration=time_duration, + sampling_frequency=sampling_frequency, duration=duration, start_time=start_time) injection_polarizations = interferometer.inject_signal( @@ -1754,7 +1754,7 @@ def get_interferometer_with_fake_noise_and_injection( def get_event_data( - event, interferometer_names=None, time_duration=4, roll_off=0.4, + event, interferometer_names=None, duration=4, roll_off=0.4, psd_offset=-1024, psd_duration=100, cache=True, outdir='outdir', label=None, plot=True, filter_freq=None, raw_data_file=None, **kwargs): """ @@ -1768,7 +1768,7 @@ def get_event_data( interferometer_names: list, optional List of interferometer identifiers, e.g., 'H1'. If None will look for data in 'H1', 'V1', 'L1' - time_duration: float + duration: float Time duration to search for. roll_off: float The roll-off (in seconds) used in the Tukey window. @@ -1808,7 +1808,7 @@ def get_event_data( for name in interferometer_names: try: interferometers.append(get_interferometer_with_open_data( - name, trigger_time=event_time, time_duration=time_duration, roll_off=roll_off, + name, trigger_time=event_time, duration=duration, roll_off=roll_off, psd_offset=psd_offset, psd_duration=psd_duration, cache=cache, outdir=outdir, label=label, plot=plot, filter_freq=filter_freq, raw_data_file=raw_data_file, **kwargs)) diff --git a/tupak/gw/likelihood.py b/tupak/gw/likelihood.py index b67ae5c449e4ac06f0ecdc82d1c07ed0faa18d07..3b197d9f17ae001a908c6bc7b53564f32d3c15c6 100644 --- a/tupak/gw/likelihood.py +++ b/tupak/gw/likelihood.py @@ -62,6 +62,7 @@ class GravitationalWaveTransient(likelihood.Likelihood): self.distance_marginalization = distance_marginalization self.phase_marginalization = phase_marginalization self.prior = prior + self._check_set_duration_and_sampling_frequency_of_waveform_generator() if self.distance_marginalization: self.check_prior_is_set() @@ -80,6 +81,27 @@ class GravitationalWaveTransient(likelihood.Likelihood): if self.time_marginalization: self.check_prior_is_set() + def _check_set_duration_and_sampling_frequency_of_waveform_generator(self): + """ Check the waveform_generator has the same duration and + sampling_frequency as the interferometers. If they are unset, then + set them, if they differ, raise an error + """ + + attributes = ['duration', 'sampling_frequency', 'start_time'] + for attr in attributes: + wfg_attr = getattr(self.waveform_generator, attr) + ifo_attr = getattr(self.interferometers, attr) + if wfg_attr is None: + logging.debug( + "The waveform_generator {} is None. Setting from the " + "provided interferometers.".format(attr)) + setattr(self.waveform_generator, attr, ifo_attr) + elif wfg_attr != ifo_attr: + logging.warning( + "The waveform_generator {} is not equal to that of the " + "provided interferometers. Overwriting the " + "waveform_generator.".format(attr)) + def check_prior_is_set(self): if self.prior is None: raise ValueError("You can't use a marginalized likelihood without specifying a prior") @@ -114,7 +136,7 @@ class GravitationalWaveTransient(likelihood.Likelihood): interferometer.frequency_domain_strain, interferometer.frequency_domain_strain, interferometer.power_spectral_density_array, - self.waveform_generator.time_duration) / 2 + self.waveform_generator.duration) / 2 return log_l.real def log_likelihood_ratio(self): @@ -130,10 +152,10 @@ class GravitationalWaveTransient(likelihood.Likelihood): signal_ifo = interferometer.get_detector_response(waveform_polarizations, self.waveform_generator.parameters) matched_filter_snr_squared += tupak.gw.utils.matched_filter_snr_squared( - signal_ifo, interferometer, self.waveform_generator.time_duration) + signal_ifo, interferometer, self.waveform_generator.duration) optimal_snr_squared += tupak.gw.utils.optimal_snr_squared( - signal_ifo, interferometer, self.waveform_generator.time_duration) + signal_ifo, interferometer, self.waveform_generator.duration) if self.time_marginalization: interferometer.time_marginalization = self.time_marginalization matched_filter_snr_squared_tc_array += 4. * (1. / interferometer.duration) * np.fft.ifft( @@ -289,7 +311,7 @@ class BasicGravitationalWaveTransient(likelihood.Likelihood): """ log_l = 0 for interferometer in self.interferometers: - log_l -= 2. / self.waveform_generator.time_duration * np.sum( + log_l -= 2. / self.waveform_generator.duration * np.sum( abs(interferometer.frequency_domain_strain) ** 2 / interferometer.power_spectral_density_array) return log_l.real @@ -330,7 +352,7 @@ class BasicGravitationalWaveTransient(likelihood.Likelihood): signal_ifo = interferometer.get_detector_response( waveform_polarizations, self.waveform_generator.parameters) - log_l = - 2. / self.waveform_generator.time_duration * np.vdot( + log_l = - 2. / self.waveform_generator.duration * np.vdot( interferometer.frequency_domain_strain - signal_ifo, (interferometer.frequency_domain_strain - signal_ifo) / interferometer.power_spectral_density_array) @@ -353,7 +375,7 @@ def get_binary_black_hole_likelihood(interferometers): """ waveform_generator = tupak.gw.waveform_generator.WaveformGenerator( - time_duration=interferometers.duration, sampling_frequency=interferometers.sampling_frequency, + 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}) return tupak.gw.likelihood.GravitationalWaveTransient(interferometers, waveform_generator) diff --git a/tupak/gw/utils.py b/tupak/gw/utils.py index fbbc5b7c16d2e7b53b9dab9a66ddbd92f71421d8..0fda7255e5108f411fd1c719559e3ca31e64b90e 100644 --- a/tupak/gw/utils.py +++ b/tupak/gw/utils.py @@ -193,7 +193,7 @@ def inner_product(aa, bb, frequency, PSD): return 4. * np.real(integral) -def noise_weighted_inner_product(aa, bb, power_spectral_density, time_duration): +def noise_weighted_inner_product(aa, bb, power_spectral_density, duration): """ Calculate the noise weighted inner product between two arrays. @@ -205,8 +205,8 @@ def noise_weighted_inner_product(aa, bb, power_spectral_density, time_duration): Array not to be complex conjugated power_spectral_density: array_like Power spectral density of the noise - time_duration: float - time_duration of the data + duration: float + duration of the data Returns ------ @@ -214,10 +214,10 @@ def noise_weighted_inner_product(aa, bb, power_spectral_density, time_duration): """ integrand = np.conj(aa) * bb / power_spectral_density - return 4 / time_duration * np.sum(integrand) + return 4 / duration * np.sum(integrand) -def matched_filter_snr_squared(signal, interferometer, time_duration): +def matched_filter_snr_squared(signal, interferometer, duration): """ Parameters @@ -226,7 +226,7 @@ def matched_filter_snr_squared(signal, interferometer, time_duration): Array containing the signal interferometer: tupak.gw.detector.Interferometer Interferometer which we want to have the data and noise from - time_duration: float + duration: float Time duration of the signal Returns @@ -236,10 +236,10 @@ def matched_filter_snr_squared(signal, interferometer, time_duration): """ return noise_weighted_inner_product( signal, interferometer.frequency_domain_strain, - interferometer.power_spectral_density_array, time_duration) + interferometer.power_spectral_density_array, duration) -def optimal_snr_squared(signal, interferometer, time_duration): +def optimal_snr_squared(signal, interferometer, duration): """ Parameters @@ -248,14 +248,14 @@ def optimal_snr_squared(signal, interferometer, time_duration): Array containing the signal interferometer: tupak.gw.detector.Interferometer Interferometer which we want to have the data and noise from - time_duration: float + duration: float Time duration of the signal Returns ------- float: The optimal signal to noise ratio possible squared """ - return noise_weighted_inner_product(signal, signal, interferometer.power_spectral_density_array, time_duration) + return noise_weighted_inner_product(signal, signal, interferometer.power_spectral_density_array, duration) def get_event_time(event): diff --git a/tupak/gw/waveform_generator.py b/tupak/gw/waveform_generator.py index f92d7392fac656cb5e15dece638761d0d7448901..fdee84a20641452a1476f98e56f56e16cfeedfab 100644 --- a/tupak/gw/waveform_generator.py +++ b/tupak/gw/waveform_generator.py @@ -6,7 +6,7 @@ import numpy as np class WaveformGenerator(object): - def __init__(self, time_duration, sampling_frequency, start_time=0, frequency_domain_source_model=None, + 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): @@ -14,9 +14,9 @@ class WaveformGenerator(object): Parameters ---------- - sampling_frequency: float + sampling_frequency: float, optional The sampling frequency - time_duration: float + duration: float, optional Time duration of data start_time: float, optional Starting time of the time array @@ -44,12 +44,12 @@ class WaveformGenerator(object): the WaveformGenerator object and initialised to `None`. """ - self.time_duration = time_duration + self.duration = duration self.sampling_frequency = sampling_frequency self.start_time = start_time self.frequency_domain_source_model = frequency_domain_source_model self.time_domain_source_model = time_domain_source_model - self.time_duration = time_duration + self.duration = duration self.sampling_frequency = sampling_frequency self.parameter_conversion = parameter_conversion self.non_standard_sampling_parameter_keys = non_standard_sampling_parameter_keys @@ -147,7 +147,7 @@ class WaveformGenerator(object): @property def frequency_array(self): - """ Frequency array for the waveforms. Automatically updates if sampling_frequency or time_duration are updated. + """ Frequency array for the waveforms. Automatically updates if sampling_frequency or duration are updated. Returns ------- @@ -156,7 +156,7 @@ class WaveformGenerator(object): if self.__frequency_array_updated is False: self.frequency_array = utils.create_frequency_series( self.sampling_frequency, - self.time_duration) + self.duration) return self.__frequency_array @frequency_array.setter @@ -166,7 +166,7 @@ class WaveformGenerator(object): @property def time_array(self): - """ Time array for the waveforms. Automatically updates if sampling_frequency or time_duration are updated. + """ Time array for the waveforms. Automatically updates if sampling_frequency or duration are updated. Returns ------- @@ -176,7 +176,7 @@ class WaveformGenerator(object): if self.__time_array_updated is False: self.__time_array = utils.create_time_series( self.sampling_frequency, - self.time_duration, + self.duration, self.start_time) self.__time_array_updated = True @@ -218,7 +218,7 @@ class WaveformGenerator(object): self.__parameters = dict.fromkeys(parameters) @property - def time_duration(self): + def duration(self): """ Allows one to set the time duration and automatically updates the frequency and time array. Returns @@ -226,11 +226,11 @@ class WaveformGenerator(object): float: The time duration. """ - return self.__time_duration + return self.__duration - @time_duration.setter - def time_duration(self, time_duration): - self.__time_duration = time_duration + @duration.setter + def duration(self, duration): + self.__duration = duration self.__frequency_array_updated = False self.__time_array_updated = False @@ -253,9 +253,9 @@ class WaveformGenerator(object): @property def start_time(self): - return self.__start_time + return self.__starting_time @start_time.setter - def start_time(self, start_time): - self.__start_time = start_time + def start_time(self, starting_time): + self.__starting_time = starting_time self.__time_array_updated = False