diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index 62317c858ed52a61b9d4432990827bc593236b05..5e304b4ecbbd71ead876598bf7231f92c2c2d86a 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -130,6 +130,7 @@ def get_polarization_tensor(ra, dec, time, psi, mode):
elif mode.lower() == 'breathing':
return np.einsum('i,j->ij', m, m) + np.einsum('i,j->ij', n, n)
+ # Calculating omega here to avoid calculation when model in [plus, cross, breathing]
omega = np.cross(m, n)
if mode.lower() == 'longitudinal':
return np.sqrt(2) * np.einsum('i,j->ij', omega, omega)
@@ -138,8 +139,7 @@ def get_polarization_tensor(ra, dec, time, psi, mode):
elif mode.lower() == 'y':
return np.einsum('i,j->ij', n, omega) + np.einsum('i,j->ij', omega, n)
else:
- logger.warning("{} not a polarization mode!".format(mode))
- return None
+ raise ValueError("{} not a polarization mode!".format(mode))
def get_vertex_position_geocentric(latitude, longitude, elevation):
@@ -380,7 +380,7 @@ def get_open_strain_data(
return strain
-def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=1, **kwargs):
+def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=0, **kwargs):
""" A function which accesses the open strain data
This uses `gwpy` to download the open data and then saves a cached copy for
@@ -416,23 +416,22 @@ def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=1
except RuntimeError:
logger.warning('Channel {} not found. Trying preset channel names'.format(channel))
- while not loaded:
- ligo_channel_types = ['GDS-CALIB_STRAIN', 'DCS-CALIB_STRAIN_C01', 'DCS-CALIB_STRAIN_C02',
- 'DCH-CLEAN_STRAIN_C02']
- virgo_channel_types = ['Hrec_hoft_V1O2Repro2A_16384Hz', 'FAKE_h_16384Hz_4R']
- channel_types = dict(H1=ligo_channel_types, L1=ligo_channel_types, V1=virgo_channel_types)
- for detector in channel_types.keys():
- for channel_type in channel_types[detector]:
- if loaded:
- break
- channel = '{}:{}'.format(detector, channel_type)
- try:
- strain = TimeSeries.read(source=file_name, channel=channel, start=start_time, end=end_time,
- **kwargs)
- loaded = True
- logger.info('Successfully read strain data for channel {}.'.format(channel))
- except RuntimeError:
- pass
+ ligo_channel_types = ['GDS-CALIB_STRAIN', 'DCS-CALIB_STRAIN_C01', 'DCS-CALIB_STRAIN_C02',
+ 'DCH-CLEAN_STRAIN_C02']
+ virgo_channel_types = ['Hrec_hoft_V1O2Repro2A_16384Hz', 'FAKE_h_16384Hz_4R']
+ channel_types = dict(H1=ligo_channel_types, L1=ligo_channel_types, V1=virgo_channel_types)
+ for detector in channel_types.keys():
+ for channel_type in channel_types[detector]:
+ if loaded:
+ break
+ channel = '{}:{}'.format(detector, channel_type)
+ try:
+ strain = TimeSeries.read(source=file_name, channel=channel, start=start_time, end=end_time,
+ **kwargs)
+ loaded = True
+ logger.info('Successfully read strain data for channel {}.'.format(channel))
+ except RuntimeError:
+ pass
if loaded:
return strain
diff --git a/test/detector_test.py b/test/detector_test.py
index edaa139fc74c12f694254986771decf16512ceee..07d107a074b9f534f1ab4bb45d747585519f4f7d 100644
--- a/test/detector_test.py
+++ b/test/detector_test.py
@@ -898,6 +898,13 @@ class TestInterferometerList(unittest.TestCase):
with self.assertRaises(ValueError):
self.ifo_list.inject_signal(injection_polarizations=None, waveform_generator=None)
+ def test_meta_data(self):
+ ifos_list = [self.ifo1, self.ifo2]
+ ifos = bilby.gw.detector.InterferometerList(ifos_list)
+ self.assertTrue(isinstance(ifos.meta_data, dict))
+ meta_data = {ifo.name: ifo.meta_data for ifo in ifos_list}
+ self.assertEqual(ifos.meta_data, meta_data)
+
@patch.object(bilby.gw.waveform_generator.WaveformGenerator, 'frequency_domain_strain')
def test_inject_signal_pol_none_calls_frequency_domain_strain(self, m):
waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
diff --git a/test/gw_utils_test.py b/test/gw_utils_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..bce6f97de0978fff982b1ee65d59e0559a8cfeab
--- /dev/null
+++ b/test/gw_utils_test.py
@@ -0,0 +1,188 @@
+from __future__ import absolute_import, division
+
+import unittest
+import os
+from shutil import rmtree
+
+import numpy as np
+import gwpy
+import lal
+import lalsimulation as lalsim
+
+import bilby
+from bilby.gw import utils as gwutils
+
+
+class TestGWUtils(unittest.TestCase):
+
+ def setUp(self):
+ self.outdir = 'outdir'
+ bilby.core.utils.check_directory_exists_and_if_not_mkdir(self.outdir)
+
+ def tearDown(self):
+ try:
+ rmtree(self.outdir)
+ except FileNotFoundError:
+ pass
+
+ def test_asd_from_freq_series(self):
+ freq_data = np.array([1, 2, 3])
+ df = 0.1
+ asd = gwutils.asd_from_freq_series(freq_data, df)
+ self.assertTrue(np.all(asd == freq_data * 2 * df ** 0.5))
+
+ def test_psd_from_freq_series(self):
+ freq_data = np.array([1, 2, 3])
+ df = 0.1
+ psd = gwutils.psd_from_freq_series(freq_data, df)
+ self.assertTrue(np.all(psd == (freq_data * 2 * df ** 0.5)**2))
+
+ def test_time_delay_from_geocenter(self):
+ det1 = np.array([0.1, 0.2, 0.3])
+ det2 = np.array([0.1, 0.2, 0.5])
+ ra = 0.5
+ dec = 0.2
+ time = 10
+ self.assertEqual(
+ gwutils.time_delay_geocentric(det1, det1, ra, dec, time), 0)
+ self.assertEqual(
+ gwutils.time_delay_geocentric(det1, det2, ra, dec, time),
+ 1.3253791114055397e-10)
+
+ def test_get_polarization_tensor(self):
+ ra = 1
+ dec = 2.0
+ time = 10
+ psi = 0.1
+ for mode in ['plus', 'cross', 'breathing', 'longitudinal', 'x', 'y']:
+ p = gwutils.get_polarization_tensor(ra, dec, time, psi, mode)
+ self.assertEqual(p.shape, (3, 3))
+ with self.assertRaises(ValueError):
+ gwutils.get_polarization_tensor(ra, dec, time, psi, 'not-a-mode')
+
+ def test_inner_product(self):
+ aa = np.array([1, 2, 3])
+ bb = np.array([5, 6, 7])
+ frequency = np.array([0.2, 0.4, 0.6])
+ PSD = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+ ip = gwutils.inner_product(aa, bb, frequency, PSD)
+ self.assertEqual(ip, 0)
+
+ def test_noise_weighted_inner_product(self):
+ aa = np.array([1e-23, 2e-23, 3e-23])
+ bb = np.array([5e-23, 6e-23, 7e-23])
+ frequency = np.array([100, 101, 102])
+ PSD = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+ psd = PSD.power_spectral_density_interpolated(frequency)
+ duration = 4
+ nwip = gwutils.noise_weighted_inner_product(aa, bb, psd, duration)
+ self.assertEqual(nwip, 239.87768033598326)
+
+ self.assertEqual(
+ gwutils.optimal_snr_squared(aa, psd, duration),
+ gwutils.noise_weighted_inner_product(aa, aa, psd, duration))
+
+ def test_matched_filter_snr(self):
+ signal = np.array([1e-23, 2e-23, 3e-23])
+ frequency_domain_strain = np.array([5e-23, 6e-23, 7e-23])
+ frequency = np.array([100, 101, 102])
+ PSD = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+ psd = PSD.power_spectral_density_interpolated(frequency)
+ duration = 4
+
+ mfsnr = gwutils.matched_filter_snr(
+ signal, frequency_domain_strain, psd, duration)
+ self.assertEqual(mfsnr, 25.510869054168282)
+
+ def test_get_event_time(self):
+ events = ['GW150914', 'GW151012', 'GW151226', 'GW170104', 'GW170608',
+ 'GW170729', 'GW170809', 'GW170814', 'GW170817', 'GW170818',
+ 'GW170823']
+ for event in events:
+ self.assertTrue(isinstance(gwutils.get_event_time(event), float))
+
+ self.assertTrue(gwutils.get_event_time('GW010290') is None)
+
+ def test_read_frame_file(self):
+ start_time = 0
+ end_time = 10
+ channel = 'H1:GDS-CALIB_STRAIN'
+ N = 100
+ times = np.linspace(start_time, end_time, N)
+ data = np.random.normal(0, 1, N)
+ ts = gwpy.timeseries.TimeSeries(data=data, times=times, t0=0)
+ ts.channel = gwpy.detector.Channel(channel)
+ ts.name = channel
+ filename = os.path.join(self.outdir, 'test.gwf')
+ ts.write(filename, format='gwf')
+
+ # Check reading without time limits
+ strain = gwutils.read_frame_file(
+ filename, start_time=None, end_time=None, channel=channel)
+ self.assertEqual(strain.channel.name, channel)
+ self.assertTrue(np.all(strain.value==data))
+
+ # Check reading with time limits
+ start_cut = 2
+ end_cut = 8
+ strain = gwutils.read_frame_file(
+ filename, start_time=start_cut, end_time=end_cut, channel=channel)
+ idxs = (times > start_cut) & (times < end_cut)
+ # Dropping the last element - for some reason gwpy drops the last element when reading in data
+ self.assertTrue(np.all(strain.value==data[idxs][:-1]))
+
+ # Check reading with unknown channels
+ strain = gwutils.read_frame_file(
+ filename, start_time=None, end_time=None)
+ self.assertTrue(np.all(strain.value==data))
+
+ # Check reading with incorrect channel
+ strain = gwutils.read_frame_file(
+ filename, start_time=None, end_time=None, channel='WRONG')
+ self.assertTrue(np.all(strain.value==data))
+
+ ts = gwpy.timeseries.TimeSeries(data=data, times=times, t0=0)
+ ts.name = 'NOT-A-KNOWN-CHANNEL'
+ ts.write(filename, format='gwf')
+ strain = gwutils.read_frame_file(
+ filename, start_time=None, end_time=None)
+ self.assertEqual(strain, None)
+
+ def test_convert_args_list_to_float(self):
+ self.assertEqual(
+ gwutils.convert_args_list_to_float(1, '2', 3.0), [1.0, 2.0, 3.0])
+ with self.assertRaises(ValueError):
+ gwutils.convert_args_list_to_float(1, '2', 'ten')
+
+ def test_lalsim_SimInspiralTransformPrecessingNewInitialConditions(self):
+ a = gwutils.lalsim_SimInspiralTransformPrecessingNewInitialConditions(
+ 0.1, 0, 0.6, 0.5, 0.6, 0.1, 0.8, 30.6, 23.2, 50, 0)
+ self.assertTrue(len(a) == 7)
+
+ def test_get_approximant(self):
+ with self.assertRaises(ValueError):
+ gwutils.lalsim_GetApproximantFromString(10)
+ self.assertEqual(gwutils.lalsim_GetApproximantFromString("IMRPhenomPV2"), 69)
+
+ def test_lalsim_SimInspiralChooseFDWaveform(self):
+ a = gwutils.lalsim_SimInspiralChooseFDWaveform(
+ 35.2, 20.4, 0.1, 0.2, 0.2, 0.2, 0.2, 0.1, 1000, 2, 2.3,
+ 45., 0.1, 10, 0.01, 10, 1000, 20, None, 69)
+ self.assertEqual(len(a), 2)
+ self.assertEqual(type(a[0]), lal.COMPLEX16FrequencySeries)
+ self.assertEqual(type(a[1]), lal.COMPLEX16FrequencySeries)
+
+ with self.assertRaises(ValueError):
+ a = gwutils.lalsim_SimInspiralChooseFDWaveform(
+ 35.2, 20.4, 0.1, 0.2, 0.2, 0.2, 0.2, 0.1, 1000, 2, 2.3,
+ 45., 0.1, 10, 0.01, 10, 1000, 20, None, 'IMRPhenomPV2')
+
+ def test_lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(self):
+ version = lalsim.IMRPhenomPv2_V
+ a = gwutils.lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(
+ 25.6, 12.2, 50, 0.2, 0, 0, 0, 0, 0, 0, 0, version)
+ self.assertEqual(len(a), 7)
+
+
+if __name__ == '__main__':
+ unittest.main()