Skip to content
Snippets Groups Projects
Commit 5eaa262b authored by Colm Talbot's avatar Colm Talbot
Browse files

Merge branch 'triangular_ifo' into 'master' 

Fix latitude and longitude calculation for triangular interferometer

See merge request !1046
parents 6004f202 fd836a86
No related branches found
No related tags found
1 merge request!1046Fix latitude and longitude calculation for triangular interferometer
Pipeline #346680 passed
Stage: initial
Stage: test
Stage: docs
Stage: deploy
Hide whitespace changes
Inline Side-by-side
bilby/gw/detector/networks.py
+ 146
− 67
View file @ 5eaa262b
......@@ -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
test/gw/detector/networks_test.py
+ 47
− 4
View file @ 5eaa262b
......@@ -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()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment