diff --git a/bilby/gw/detector/networks.py b/bilby/gw/detector/networks.py index 4370219c183e2a20e86ad53fece558cdd147f3a4..bc1067ff94b12ff0a15bccc6c7773f3fa6c6cd37 100644 --- a/bilby/gw/detector/networks.py +++ b/bilby/gw/detector/networks.py @@ -11,10 +11,10 @@ from .psd import PowerSpectralDensity class InterferometerList(list): - """ A list of Interferometer objects """ + """A list of Interferometer objects""" def __init__(self, interferometers): - """ Instantiate a InterferometerList + """Instantiate a InterferometerList The InterferometerList is a list of Interferometer objects, each object has the data used in evaluating the likelihood @@ -32,7 +32,9 @@ class InterferometerList(list): if type(ifo) == str: ifo = get_empty_interferometer(ifo) if type(ifo) not in [Interferometer, TriangularInterferometer]: - raise TypeError("Input list of interferometers are not all Interferometer objects") + raise TypeError( + "Input list of interferometers are not all Interferometer objects" + ) else: self.append(ifo) self._check_interferometers() @@ -45,28 +47,37 @@ class InterferometerList(list): If both checks fail, then a ValueError is raised. """ - consistent_attributes = ['duration', 'start_time', 'sampling_frequency'] + consistent_attributes = ["duration", "start_time", "sampling_frequency"] for attribute in consistent_attributes: - x = [getattr(interferometer.strain_data, attribute) - for interferometer in self] + x = [ + getattr(interferometer.strain_data, attribute) + for interferometer in self + ] try: if not all(y == x[0] for y in x): - ifo_strs = ["{ifo}[{attribute}]={value}".format( - ifo=ifo.name, - attribute=attribute, - value=getattr(ifo.strain_data, attribute)) - for ifo in self] + ifo_strs = [ + "{ifo}[{attribute}]={value}".format( + ifo=ifo.name, + attribute=attribute, + value=getattr(ifo.strain_data, attribute), + ) + for ifo in self + ] raise ValueError( "The {} of all interferometers are not the same: {}".format( - attribute, ', '.join(ifo_strs))) + attribute, ", ".join(ifo_strs) + ) + ) except ValueError as e: if not all(math.isclose(y, x[0], abs_tol=1e-5) for y in x): raise ValueError(e) else: logger.warning(e) - def set_strain_data_from_power_spectral_densities(self, sampling_frequency, duration, start_time=0): - """ Set the `Interferometer.strain_data` from the power spectral densities of the detectors + def set_strain_data_from_power_spectral_densities( + self, sampling_frequency, duration, start_time=0 + ): + """Set the `Interferometer.strain_data` from the power spectral densities of the detectors This uses the `interferometer.power_spectral_density` object to set the `strain_data` to a noise realization. See @@ -83,12 +94,16 @@ class InterferometerList(list): """ for interferometer in self: - interferometer.set_strain_data_from_power_spectral_density(sampling_frequency=sampling_frequency, - duration=duration, - start_time=start_time) + interferometer.set_strain_data_from_power_spectral_density( + sampling_frequency=sampling_frequency, + duration=duration, + start_time=start_time, + ) - def set_strain_data_from_zero_noise(self, sampling_frequency, duration, start_time=0): - """ Set the `Interferometer.strain_data` from the power spectral densities of the detectors + def set_strain_data_from_zero_noise( + self, sampling_frequency, duration, start_time=0 + ): + """Set the `Interferometer.strain_data` from the power spectral densities of the detectors This uses the `interferometer.power_spectral_density` object to set the `strain_data` to zero noise. See @@ -105,9 +120,11 @@ class InterferometerList(list): """ for interferometer in self: - interferometer.set_strain_data_from_zero_noise(sampling_frequency=sampling_frequency, - duration=duration, - start_time=start_time) + interferometer.set_strain_data_from_zero_noise( + sampling_frequency=sampling_frequency, + duration=duration, + start_time=start_time, + ) def inject_signal( self, @@ -147,12 +164,14 @@ class InterferometerList(list): """ if injection_polarizations is None: if waveform_generator is not None: - injection_polarizations = \ - waveform_generator.frequency_domain_strain(parameters) + injection_polarizations = waveform_generator.frequency_domain_strain( + parameters + ) else: raise ValueError( "inject_signal needs one of waveform_generator or " - "injection_polarizations.") + "injection_polarizations." + ) all_injection_polarizations = list() for interferometer in self: @@ -167,7 +186,7 @@ class InterferometerList(list): return all_injection_polarizations def save_data(self, outdir, label=None): - """ Creates a save file for the data in plain text format + """Creates a save file for the data in plain text format Parameters ========== @@ -179,7 +198,7 @@ class InterferometerList(list): for interferometer in self: interferometer.save_data(outdir=outdir, label=label) - def plot_data(self, signal=None, outdir='.', label=None): + def plot_data(self, signal=None, outdir=".", label=None): if utils.command_line_args.bilby_test_mode: return @@ -223,13 +242,14 @@ class InterferometerList(list): @property def meta_data(self): - """ Dictionary of the per-interferometer meta_data """ - return {interferometer.name: interferometer.meta_data - for interferometer in self} + """Dictionary of the per-interferometer meta_data""" + return { + interferometer.name: interferometer.meta_data for interferometer in self + } @staticmethod def _filename_from_outdir_label_extension(outdir, label, extension="h5"): - return os.path.join(outdir, label + f'.{extension}') + return os.path.join(outdir, label + f".{extension}") _save_docstring = """ Saves the object to a {format} file @@ -253,53 +273,66 @@ class InterferometerList(list): """ - def to_hdf5(self, outdir='outdir', label='ifo_list'): + def to_hdf5(self, outdir="outdir", label="ifo_list"): import deepdish + if sys.version_info[0] < 3: - raise NotImplementedError('Pickling of InterferometerList is not supported in Python 2.' - 'Use Python 3 instead.') - label = label + '_' + ''.join(ifo.name for ifo in self) + raise NotImplementedError( + "Pickling of InterferometerList is not supported in Python 2." + "Use Python 3 instead." + ) + label = label + "_" + "".join(ifo.name for ifo in self) utils.check_directory_exists_and_if_not_mkdir(outdir) try: filename = self._filename_from_outdir_label_extension(outdir, label, "h5") deepdish.io.save(filename, self) except AttributeError: - logger.warning("Saving to hdf5 using deepdish failed. Pickle dumping instead.") + logger.warning( + "Saving to hdf5 using deepdish failed. Pickle dumping instead." + ) self.to_pickle(outdir=outdir, label=label) @classmethod def from_hdf5(cls, filename=None): import deepdish + if sys.version_info[0] < 3: - raise NotImplementedError('Pickling of InterferometerList is not supported in Python 2.' - 'Use Python 3 instead.') + raise NotImplementedError( + "Pickling of InterferometerList is not supported in Python 2." + "Use Python 3 instead." + ) res = deepdish.io.load(filename) if res.__class__ == list: res = cls(res) if res.__class__ != cls: - raise TypeError('The loaded object is not a InterferometerList') + raise TypeError("The loaded object is not a InterferometerList") return res def to_pickle(self, outdir="outdir", label="ifo_list"): import dill - utils.check_directory_exists_and_if_not_mkdir('outdir') - label = label + '_' + ''.join(ifo.name for ifo in self) - filename = self._filename_from_outdir_label_extension(outdir, label, extension="pkl") + + utils.check_directory_exists_and_if_not_mkdir("outdir") + label = label + "_" + "".join(ifo.name for ifo in self) + filename = self._filename_from_outdir_label_extension( + outdir, label, extension="pkl" + ) with open(filename, "wb") as ff: dill.dump(self, ff) @classmethod def from_pickle(cls, filename=None): import dill + with open(filename, "rb") as ff: res = dill.load(ff) if res.__class__ != cls: - raise TypeError('The loaded object is not an InterferometerList') + raise TypeError("The loaded object is not an InterferometerList") return res to_hdf5.__doc__ = _save_docstring.format( - format="hdf5", extra=""".. deprecated:: 1.1.0 - Use :func:`to_pickle` instead.""" + format="hdf5", + extra=""".. deprecated:: 1.1.0 + Use :func:`to_pickle` instead.""", ) to_pickle.__doc__ = _save_docstring.format( format="pickle", extra=".. versionadded:: 1.1.0" @@ -309,10 +342,21 @@ class InterferometerList(list): class TriangularInterferometer(InterferometerList): - - def __init__(self, name, power_spectral_density, minimum_frequency, maximum_frequency, - length, latitude, longitude, elevation, xarm_azimuth, yarm_azimuth, - xarm_tilt=0., yarm_tilt=0.): + def __init__( + self, + name, + power_spectral_density, + minimum_frequency, + maximum_frequency, + length, + latitude, + longitude, + elevation, + xarm_azimuth, + yarm_azimuth, + xarm_tilt=0.0, + yarm_tilt=0.0, + ): super(TriangularInterferometer, self).__init__([]) self.name = name # for attr in ['power_spectral_density', 'minimum_frequency', 'maximum_frequency']: @@ -324,15 +368,46 @@ class TriangularInterferometer(InterferometerList): maximum_frequency = [maximum_frequency] * 3 for ii in range(3): - self.append(Interferometer( - '{}{}'.format(name, ii + 1), power_spectral_density[ii], minimum_frequency[ii], maximum_frequency[ii], - length, latitude, longitude, elevation, xarm_azimuth, yarm_azimuth, xarm_tilt, yarm_tilt)) + self.append( + Interferometer( + "{}{}".format(name, ii + 1), + power_spectral_density[ii], + minimum_frequency[ii], + maximum_frequency[ii], + length, + latitude, + longitude, + elevation, + xarm_azimuth, + yarm_azimuth, + xarm_tilt, + yarm_tilt, + ) + ) xarm_azimuth += 240 yarm_azimuth += 240 - latitude += np.arctan(length * np.sin(xarm_azimuth * np.pi / 180) * 1e3 / utils.radius_of_earth) - longitude += np.arctan(length * np.cos(xarm_azimuth * np.pi / 180) * 1e3 / utils.radius_of_earth) + latitude += ( + np.arctan( + length + * np.sin(xarm_azimuth * np.pi / 180) + * 1e3 + / utils.radius_of_earth + ) + * 180 + / np.pi + ) + longitude += ( + np.arctan( + length + * np.cos(xarm_azimuth * np.pi / 180) + * 1e3 + / utils.radius_of_earth + ) + * 180 + / np.pi + ) def get_empty_interferometer(name): @@ -365,34 +440,38 @@ def get_empty_interferometer(name): interferometer: Interferometer Interferometer instance """ - filename = os.path.join(os.path.dirname(__file__), 'detectors', '{}.interferometer'.format(name)) + filename = os.path.join( + os.path.dirname(__file__), "detectors", "{}.interferometer".format(name) + ) try: return load_interferometer(filename) except OSError: - raise ValueError('Interferometer {} not implemented'.format(name)) + raise ValueError("Interferometer {} not implemented".format(name)) def load_interferometer(filename): """Load an interferometer from a file.""" parameters = dict() - with open(filename, 'r') as parameter_file: + with open(filename, "r") as parameter_file: lines = parameter_file.readlines() for line in lines: - if line[0] == '#' or line[0] == '\n': + if line[0] == "#" or line[0] == "\n": continue - split_line = line.split('=') + split_line = line.split("=") key = split_line[0].strip() - value = eval('='.join(split_line[1:])) + value = eval("=".join(split_line[1:])) parameters[key] = value - if 'shape' not in parameters.keys(): + if "shape" not in parameters.keys(): ifo = Interferometer(**parameters) - logger.debug('Assuming L shape for {}'.format('name')) - elif parameters['shape'].lower() in ['l', 'ligo']: - parameters.pop('shape') + logger.debug("Assuming L shape for {}".format("name")) + elif parameters["shape"].lower() in ["l", "ligo"]: + parameters.pop("shape") ifo = Interferometer(**parameters) - elif parameters['shape'].lower() in ['triangular', 'triangle']: - parameters.pop('shape') + elif parameters["shape"].lower() in ["triangular", "triangle"]: + parameters.pop("shape") ifo = TriangularInterferometer(**parameters) else: - raise IOError("{} could not be loaded. Invalid parameter 'shape'.".format(filename)) + raise IOError( + "{} could not be loaded. Invalid parameter 'shape'.".format(filename) + ) return ifo diff --git a/test/gw/detector/networks_test.py b/test/gw/detector/networks_test.py index 2ad4060e5934a571e9b902eb22119a9303a80cf9..b0d02e4eae2e58f7ad43aee61be3365c86137603 100644 --- a/test/gw/detector/networks_test.py +++ b/test/gw/detector/networks_test.py @@ -3,6 +3,7 @@ import unittest import pytest from shutil import rmtree from packaging import version +from itertools import combinations import deepdish as dd import mock @@ -119,7 +120,7 @@ class TestInterferometerList(unittest.TestCase): bilby.gw.detector.InterferometerList("string") def test_init_with_string_list(self): - """ Merely checks if this ends up in the right bracket """ + """Merely checks if this ends up in the right bracket""" with mock.patch("bilby.gw.detector.networks.get_empty_interferometer") as m: m.side_effect = TypeError with self.assertRaises(TypeError): @@ -240,8 +241,12 @@ class TestInterferometerList(unittest.TestCase): @patch.object(bilby.gw.detector.Interferometer, "inject_signal") def test_inject_signal_with_inj_pol(self, m): - self.ifo_list.inject_signal(injection_polarizations=dict(plus=1), raise_error=False) - m.assert_called_with(parameters=None, injection_polarizations=dict(plus=1), raise_error=False) + self.ifo_list.inject_signal( + injection_polarizations=dict(plus=1), raise_error=False + ) + m.assert_called_with( + parameters=None, injection_polarizations=dict(plus=1), raise_error=False + ) self.assertEqual(len(self.ifo_list), m.call_count) @patch.object(bilby.gw.detector.Interferometer, "inject_signal") @@ -337,7 +342,9 @@ class TestInterferometerList(unittest.TestCase): recovered_ifo = bilby.gw.detector.InterferometerList.from_hdf5(filename) self.assertListEqual(self.ifo_list, recovered_ifo) - @pytest.mark.skipif(pandas_version_test or sys.version_info[0] < 3, reason=skip_reason) + @pytest.mark.skipif( + pandas_version_test or sys.version_info[0] < 3, reason=skip_reason + ) def test_to_and_from_hdf5_wrong_class(self): dd.io.save("./outdir/psd.h5", self.ifo_list[0].power_spectral_density) filename = self.ifo_list._filename_from_outdir_label_extension( @@ -354,6 +361,7 @@ class TestInterferometerList(unittest.TestCase): def test_to_and_from_pkl_wrong_class(self): import dill + with open("./outdir/psd.pkl", "wb") as ff: dill.dump(self.ifo_list[0].power_spectral_density, ff) filename = self.ifo_list._filename_from_outdir_label_extension( @@ -372,5 +380,40 @@ class TestInterferometerList(unittest.TestCase): ifos.plot_data(outdir=self.outdir) +class TriangularInterferometerTest(unittest.TestCase): + def setUp(self): + self.triangular_ifo = bilby.gw.detector.get_empty_interferometer("ET") + + def tearDown(self): + del self.triangular_ifo + + def test_individual_positions(self): + """ + Check that the distances between the vertices of the three + individual interferometers is approximately equal to the + length of the arms of the detector. + Calculation following: + https://www.movable-type.co.uk/scripts/latlong.html + Angles must be in radians + """ + + def a(delta_lat, delta_long, lat_1, lat_2): + return ( + np.sin(delta_lat / 2) ** 2 + + np.cos(lat_1) * np.cos(lat_2) * np.sin(delta_long / 2) ** 2 + ) + + def c(a): + return 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a)) + + for pair in list(combinations(self.triangular_ifo, 2)): + delta_lat = np.radians(pair[1].latitude - pair[0].latitude) + delta_long = np.radians(pair[1].longitude - pair[0].longitude) + pair_a = a(delta_lat, delta_long, pair[0].latitude, pair[1].latitude) + pair_c = c(pair_a) + distance = bilby.core.utils.radius_of_earth * pair_c + self.assertAlmostEqual(distance / 1000, pair[0].length, delta=1) + + if __name__ == "__main__": unittest.main()