An inject_glitch method for the interferometer object

README.rst

@@ -1,44 +1,21 @@
-|pipeline status| |coverage report| |pypi| |conda| |version|
+
 
 =====
 Bilby
 =====
 
-A user-friendly Bayesian inference library.
-Fulfilling all your Bayesian dreams.
-
-Online material to help you get started:
-
--  `Installation instructions <https://bilby-dev.github.io/bilby/installation.html>`__
--  `Documentation <https://bilby-dev.github.io/bilby/>`__
-
-If you need help, find an issue, or just have a question/suggestion you can
-
-- Join our `Slack workspace <https://bilby-code.slack.com/>`__ via this `invite link <https://join.slack.com/t/bilby-code/shared_invite/zt-2s5a0jy1g-xB7uIy1fGGxW0CkBwwTbWQp>`__
-- Ask questions (or search through other users questions and answers) on `StackOverflow <https://stackoverflow.com/questions/tagged/bilby>`__ using the bilby tag
-- Submit issues directly through `the issue tracker <https://github.com/bilby-dev/bilby/issues>`__
-- For chat.ligo.org users, join the `#bilby-help <https://chat.ligo.org/ligo/channels/bilby-help>`__ or `#bilby-devel <https://chat.ligo.org/ligo/channels/bilby-devel>`__ channels
-- For LVK-confidential issues, please open `a confidential issue on bilby_pipe <https://git.ligo.org/lscsoft/bilby_pipe/-/issues/new>`__
-
-We encourage you to contribute to the development of bilby. This is done via pull request.  For
-help in creating a pull request, see `this page
-<https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/proposing-changes-to-your-work-with-pull-requests/creating-a-pull-request>`__ or contact
-us directly. For advice on contributing, see `the contributing guide <https://github.com/bilby-dev/bilby/blob/main/CONTRIBUTING.md>`__.
-
-
---------------
-Citation guide
---------------
-
-Please refer to the `Acknowledging/citing bilby guide <https://bilby-dev.github.io/bilby/citing-bilby.html>`__.
-
-.. |pipeline status| image:: https://github.com/bilby-dev/bilby/actions/workflows/unit-tests.yml/badge.svg
-   :target: https://github.com/bilby-dev/bilby/commits/master
-.. |coverage report| image:: https://codecov.io/gh/bilby-dev/bilby/graph/badge.svg?token=BpfcD0hFSu
-   :target: https://codecov.io/gh/bilby-dev/bilby
-.. |pypi| image:: https://badge.fury.io/py/bilby.svg
-   :target: https://pypi.org/project/bilby/
-.. |conda| image:: https://img.shields.io/conda/vn/conda-forge/bilby.svg
-   :target: https://anaconda.org/conda-forge/bilby
-.. |version| image:: https://img.shields.io/pypi/pyversions/bilby.svg
-   :target: https://pypi.org/project/bilby/
+A bilby fork to develope features relevant for A# and/or 3G detectors:
+
+- Inject glitches to the data:
+   - An example script is at ``examples/gw_examples/injection_examples/glitch.py``
+   - An example script to inject `blip glitches <https://git.ligo.org/melissa.lopez/gengli>`_ is at: ``examples/gw_examples/injection_examples/inject_glitch_from_gengli.py``. 
+   - This feature is implemented by adding an ``inject_glitch`` method to the ``interferometer.py``. 
+   - Useful to estimate the impact of glitches on PE or glitch mitigation related simulations.
+
+- Earth rotation:
+   - An example script is at ``examples/gw_examples/injection_examples/rotation.py``
+   - This feature is implemented by adding a bool (``earth_rotation``) argument to the ``inject_signal`` method.
+   - Also added the same bool to the base likelihood object, to control whether both (injection/recovery) should account for Earth rotation.
+   - Currently only supports ``geocent_time`` as the reference time.
+   - Does not support time marginalisation.
+

bilby/gw/detector/networks.py

@@ -140,6 +140,7 @@ def inject_signal(
         injection_polarizations=None,
         waveform_generator=None,
         raise_error=True,
+        earth_rotation=False
     ):
         """ Inject a signal into noise in each of the three detectors.
 
@@ -188,6 +189,7 @@ def inject_signal(
                     parameters=parameters,
                     injection_polarizations=injection_polarizations,
                     raise_error=raise_error,
+                    earth_rotation=earth_rotation
                 )
             )
 

bilby/gw/likelihood/base.py

@@ -147,7 +147,7 @@ def __init__(
             distance_marginalization=False, phase_marginalization=False, calibration_marginalization=False, priors=None,
             distance_marginalization_lookup_table=None, calibration_lookup_table=None,
             number_of_response_curves=1000, starting_index=0, jitter_time=True, reference_frame="sky",
-            time_reference="geocenter"
+            time_reference="geocenter", earth_rotation=False
     ):
 
         self.waveform_generator = waveform_generator
@@ -162,6 +162,13 @@ def __init__(
         self._noise_log_likelihood_value = None
         self.jitter_time = jitter_time
         self.reference_frame = reference_frame
+        self.earth_rotation = earth_rotation
+
+        if self.earth_rotation and self.time_marginalization:
+            raise ValueError(
+                "earth_rotation and time_marginalization cannot be used "
+                "together.")
+
         if "geocent" not in time_reference:
             self.time_reference = time_reference
             self.reference_ifo = get_empty_interferometer(self.time_reference)
@@ -734,7 +741,8 @@ def _compute_full_waveform(self, signal_polarizations, interferometer, parameter
             Interferometer to compute the response with respect to.
         """
         parameters = _fallback_to_parameters(self, parameters)
-        return interferometer.get_detector_response(signal_polarizations, parameters)
+        return interferometer.get_detector_response(signal_polarizations,
+                                                    parameters, earth_rotation=self.earth_rotation)
 
     def generate_phase_sample_from_marginalized_likelihood(
             self, signal_polarizations=None, parameters=None):

examples/gw_examples/injection_examples/glitch.py

@@ -0,0 +1,139 @@
+import bilby
+import matplotlib.pyplot as plt
+import numpy as np
+
+duration = 8.0
+sampling_frequency = 2048.0
+glitch_snr = 12
+injection_parameters = dict(
+    mass_1=36.0,
+    mass_2=29.0,
+    a_1=0.4,
+    a_2=0.3,
+    tilt_1=0.5,
+    tilt_2=1.0,
+    phi_12=1.7,
+    phi_jl=0.3,
+    luminosity_distance=2100.0,
+    theta_jn=0.4,
+    psi=2.659,
+    phase=1.3,
+    geocent_time=205.0,
+    ra=1.375,
+    dec=-1.2108,
+)
+
+bilby.core.utils.random.seed(23)
+
+
+def sine_gaussian_glitch(
+    time_array, amplitude=1e-22, f0=100.0, tau=0.1, t0=0.0, phase=0.0
+):
+    """Generate a sine-Gaussian waveform.
+
+    Parameters
+    ----------
+    time_array : array_like
+        Time samples.
+    amplitude : float
+        Peak amplitude.
+    f0 : float
+        Central frequency in Hz.
+    tau : float
+        Decay time (Gaussian width) in seconds.
+    t0 : float
+        Central time of the burst in seconds.
+    phase : float
+        Phase offset in radians.
+    """
+    envelope = np.exp(-((time_array - t0) ** 2) / (2 * tau**2))
+    return amplitude * envelope * np.sin(2 * np.pi * f0 * (time_array - t0) + phase)
+
+
+# Fixed arguments passed into the source model
+waveform_arguments = dict(
+    waveform_approximant="IMRPhenomXPHM",
+    reference_frequency=50.0,
+    minimum_frequency=20.0,
+)
+
+# Create the waveform_generator using a LAL BinaryBlackHole source function
+# the generator will convert all the parameters
+waveform_generator = bilby.gw.WaveformGenerator(
+    duration=duration,
+    sampling_frequency=sampling_frequency,
+    frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+    parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+    waveform_arguments=waveform_arguments,
+)
+
+# Set up interferometers.  In this case we'll use two interferometers
+# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design
+# sensitivity
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1"])
+ifos.set_strain_data_from_zero_noise(
+    sampling_frequency=sampling_frequency,
+    duration=duration,
+    start_time=injection_parameters["geocent_time"] - 4,
+)
+
+ifos.inject_signal(
+    waveform_generator=waveform_generator, parameters=injection_parameters
+)
+
+inject_glitch_from_time_domain_strain = True
+
+if inject_glitch_from_time_domain_strain:
+    glitch_sample_times = np.arange(0, 2, 1 / 256)
+    from scipy.signal.windows import tukey
+
+    glitch = sine_gaussian_glitch(glitch_sample_times, tau=0.1, t0=0.5) * tukey(
+        len(glitch_sample_times), alpha=0.2
+    )
+
+    # Inject glitch
+    glitch_parameters = {
+        "onset_time": injection_parameters["geocent_time"] + 1,
+        "snr": glitch_snr,
+    }
+    ifos[0].inject_glitch(
+        glitch_snr,
+        glitch_parameters=glitch_parameters,
+        glitch_time_domain_strain=glitch,
+        glitch_sample_times=glitch_sample_times,
+    )
+
+
+inject_glitch_from_a_glitch_waveform_generator = False
+if inject_glitch_from_a_glitch_waveform_generator:
+
+    glitch_waveform_generator = bilby.gw.WaveformGenerator(
+        duration=duration,
+        sampling_frequency=sampling_frequency,
+        time_domain_source_model=sine_gaussian_glitch,
+    )
+    glitch_parameters = {
+        "onset_time": injection_parameters["geocent_time"] + 1,
+        "amplitude": 1e-22,
+        "f0": 100.0,
+        "tau": 0.1,
+        "t0": 1.3,
+        "phase": 0.0,
+        "snr": glitch_snr,
+    }
+
+    ifos[0].inject_glitch(
+        glitch_parameters=glitch_parameters,
+        glitch_waveform_generator=glitch_waveform_generator,
+    )
+
+
+fig, axes = plt.subplots(2, 1)
+ax = axes[0]
+
+ax.axvline(x=injection_parameters["geocent_time"], color="black", label="Merger time")
+ax.axvline(x=glitch_parameters["onset_time"], color="blue", label="glitch onset time")
+ax.legend()
+
+ax.plot(ifos[0].time_array, ifos[0].time_domain_strain)
+fig.savefig("./time_domain_strain.pdf")

examples/gw_examples/injection_examples/inject_glitch_from_gengli.py

@@ -0,0 +1,88 @@
+import bilby
+import gengli
+import matplotlib.pyplot as plt
+import numpy as np
+
+# gengli package needs to be installed from ``pip install gengli``
+
+duration = 8.0
+sampling_frequency = 2048.0
+glitch_snr = 12
+injection_parameters = dict(
+    mass_1=36.0,
+    mass_2=29.0,
+    a_1=0.4,
+    a_2=0.3,
+    tilt_1=0.5,
+    tilt_2=1.0,
+    phi_12=1.7,
+    phi_jl=0.3,
+    luminosity_distance=2100.0,
+    theta_jn=0.4,
+    psi=2.659,
+    phase=1.3,
+    geocent_time=205.0,
+    ra=1.375,
+    dec=-1.2108,
+)
+
+bilby.core.utils.random.seed(23)
+
+# Fixed arguments passed into the source model
+waveform_arguments = dict(
+    waveform_approximant="IMRPhenomXPHM",
+    reference_frequency=50.0,
+    minimum_frequency=20.0,
+)
+
+# Create the waveform_generator using a LAL BinaryBlackHole source function
+# the generator will convert all the parameters
+waveform_generator = bilby.gw.WaveformGenerator(
+    duration=duration,
+    sampling_frequency=sampling_frequency,
+    frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+    parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+    waveform_arguments=waveform_arguments,
+)
+
+# Set up interferometers.  In this case we'll use two interferometers
+# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design
+# sensitivity
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1"])
+ifos.set_strain_data_from_zero_noise(
+    sampling_frequency=sampling_frequency,
+    duration=duration,
+    start_time=injection_parameters["geocent_time"] - 4,
+)
+
+ifos.inject_signal(
+    waveform_generator=waveform_generator, parameters=injection_parameters
+)
+
+glitch_generator = gengli.glitch_generator("L1")
+
+# Produces whitened time domain blip glitch
+glitch_time_domain_strain = glitch_generator.get_glitch(srate=sampling_frequency)
+
+# onset_time is the time at which glitch will be injected and snr is the snr of the glitch.
+glitch_parameters = dict(onset_time=injection_parameters["geocent_time"] + 2, snr=9)
+glitch_sample_times = (
+    np.arange(0, len(glitch_time_domain_strain), 1) / sampling_frequency
+)
+
+ifos[0].inject_glitch(
+    glitch_snr,
+    glitch_parameters=glitch_parameters,
+    glitch_time_domain_strain=glitch_time_domain_strain,
+    glitch_sample_times=glitch_sample_times,
+)
+
+fig, axes = plt.subplots(2, 1)
+ax = axes[0]
+
+ax.axvline(x=injection_parameters["geocent_time"], color="black", label="Merger time")
+ax.axvline(x=glitch_parameters["onset_time"], color="blue", label="glitch onset time")
+ax.legend()
+
+ax.plot(ifos[0].time_array, ifos[0].time_domain_strain)
+fig.savefig("./time_domain_strain.pdf")