Skip to content
Snippets Groups Projects
Commit 347418ac authored by Moritz Huebner's avatar Moritz Huebner
Browse files

Merge branch 'master' into tests_update

# Conflicts:
#	tupak/gw/waveform_generator.py
parents 55229b01 2ca87612
No related branches found
No related tags found
No related merge requests found
Showing
with 132 additions and 95 deletions
......@@ -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
......
......@@ -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
......
......@@ -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']]
......
......@@ -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']]
......
......@@ -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
......
......@@ -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
......
......@@ -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
......
......@@ -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})
......
......@@ -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)
......
......@@ -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):
......
......@@ -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]
......
......@@ -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))
......
......@@ -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)
......
......@@ -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):
......
......@@ -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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment