FEAT: add argument to specify calibration correction direction (#47)

* FEAT: add argument to specify calibration correction direction

* FEAT: add calibration correction specification to the calibration draws file

* FORMAT: flake fix in calibration

* BUG: fix calibration table read/write

* DOC: Improve documentation for calibration correction types

* DOC: add module-level docstring describing calibration conventions

* TYPO: refactoring typo fixes

* TYPO: typo fix in documentation

* TEST: update argument name in test

* DOC: fix versionadded

* DOC: fix version changed

bilby/gw/detector/calibration.py

@@ -1,4 +1,43 @@
-""" Functions for adding calibration factors to waveform templates.
+r"""
+Functions for adding calibration factors to waveform templates.
+
+The two key quantities are :math:`d`, the (possible mis-)calibrated strain
+data used for the analysis, and :math:`h`, the theoretical strain predicted by
+the waveform model.
+
+There are two conventions in the literature for how to specify calibration
+corrections. People who work on gravitational-wave detector calibration
+typically describe a correction to the data so that the signal matches the
+theoretical prediction
+
+.. math::
+
+    h = \eta d.
+
+However, when performing inference, we are interested in the correction that
+must be applied to the theoretical strain to match the signal contained within
+the data
+
+.. math::
+
+    d = \alpha h.
+
+Clearly, these are related via
+
+.. math::
+
+    \eta = \frac{1}{\alpha}
+
+Internally, in :code:`Bilby`, the correction is always :math:`\alpha`.
+However, when reading in a data product describing calibration uncertainty,
+e.g., uncertainty envelopes or estimated response curves, the user should
+specify which method is being used as :code:`"data"` for :math:`\eta` or
+:code:`"template"` for :math:`\alpha`.
+
+.. note::
+    In general, data products produced by the LVK calibration groups use the
+    :code:`"data"` convention.
+
 """
 import copy
 import os
@@ -12,9 +51,34 @@
 from ..prior import CalibrationPriorDict
 
 
-def read_calibration_file(filename, frequency_array, number_of_response_curves, starting_index=0):
-    """
-    Function to read the hdf5 files from the calibration group containing the physical calibration response curves.
+def _check_calibration_correction_type(correction_type):
+    if correction_type is None:
+        logger.warning(
+            "Calibration envelope correction type is not specified. "
+            "Assuming this correction maps calibrated data to theoretical "
+            "strain. If this is correct, this should be explicitly "
+            "specified via CalibrationPriorDict.from_envelope_file(..., "
+            "correction_type='data'). The other possibility is correction_type="
+            "'template', which maps theoretical strain to calibrated data."
+        )
+        correction_type = "data"
+    if correction_type.lower() not in ["data", "template"]:
+        raise ValueError(
+            "Calibration envelope correction should be one of 'data' or "
+            f"'template', found {correction_type}."
+        )
+    logger.debug(
+        f"Supplied calibration correction will be applied to the {correction_type}"
+    )
+    return correction_type
+
+
+def read_calibration_file(
+    filename, frequency_array, number_of_response_curves, starting_index=0, correction_type=None
+):
+    r"""
+    Function to read the hdf5 files from the calibration group containing the
+    physical calibration response curves.
 
     Parameters
     ----------
@@ -27,6 +91,14 @@ def read_calibration_file(filename, frequency_array, number_of_response_curves,
     starting_index: int
         Index of the first curve to use within the array. This allows for segmenting the calibration curve array
         into smaller pieces.
+    correction_type: str
+        How the correction is defined, either to the :code:`data`
+        (default) or the :code:`template`. In general, data products
+        produced by the LVK calibration groups assume :code:`data`.
+        The default value will be removed in a future release and
+        this will need to be explicitly specified.
+
+        .. versionadded:: 1.4.0
 
     Returns
     -------
@@ -37,6 +109,8 @@ def read_calibration_file(filename, frequency_array, number_of_response_curves,
     """
     import tables
 
+    correction_type = _check_calibration_correction_type(correction_type=correction_type)
+
     logger.info(f"Reading calibration draws from {filename}")
     calibration_file = tables.open_file(filename, 'r')
     calibration_amplitude = \
@@ -58,6 +132,8 @@ def read_calibration_file(filename, frequency_array, number_of_response_curves,
     calibration_draws = interp1d(
         calibration_frequencies, calibration_draws, kind='cubic',
         bounds_error=False, fill_value=1)(frequency_array)
+    if correction_type == "data":
+        calibration_draws = 1 / calibration_draws
 
     try:
         parameter_draws = pd.read_hdf(filename, key="CalParams")
@@ -67,7 +143,9 @@ def read_calibration_file(filename, frequency_array, number_of_response_curves,
     return calibration_draws, parameter_draws
 
 
-def write_calibration_file(filename, frequency_array, calibration_draws, calibration_parameter_draws=None):
+def write_calibration_file(
+    filename, frequency_array, calibration_draws, calibration_parameter_draws=None, correction_type=None
+):
     """
     Function to write the generated response curves to file
 
@@ -82,10 +160,23 @@ def write_calibration_file(filename, frequency_array, calibration_draws, calibra
         Shape is (number_of_response_curves x len(frequency_array))
     calibration_parameter_draws: data_frame
         Parameters used to generate the random draws of the calibration response curves
+    correction_type: str
+        How the correction is defined, either to the :code:`data`
+        (default) or the :code:`template`. In general, data products
+        produced by the LVK calibration groups assume :code:`data`.
+        The default value will be removed in a future release and
+        this will need to be explicitly specified.
+
+        .. versionadded:: 1.4.0
 
     """
     import tables
 
+    correction_type = _check_calibration_correction_type(correction_type=correction_type)
+
+    if correction_type == "data":
+        calibration_draws = 1 / calibration_draws
+
     logger.info(f"Writing calibration draws to {filename}")
     calibration_file = tables.open_file(filename, 'w')
     deltaR_group = calibration_file.create_group(calibration_file.root, 'deltaR')

bilby/gw/prior.py

@@ -1138,7 +1138,7 @@ def to_file(self, outdir, label):
     @staticmethod
     def from_envelope_file(envelope_file, minimum_frequency,
                            maximum_frequency, n_nodes, label,
-                           boundary="reflective"):
+                           boundary="reflective", correction_type=None):
         """
         Load in the calibration envelope.
 
@@ -1149,6 +1149,14 @@ def from_envelope_file(envelope_file, minimum_frequency,
             frequency median-amplitude median-phase -1-sigma-amplitude
             -1-sigma-phase +1-sigma-amplitude +1-sigma-phase
 
+        There are two definitions of the calibration correction in the
+        literature, one defines the correction as mapping calibrated strain
+        to theoretical waveform templates (:code:`data`) and the other as
+        mapping theoretical waveform templates to calibrated strain
+        (:code:`template`). Prior to version 1.4.0, :code:`template` was assumed,
+        the default changed to :code:`data` when the :code:`correction` argument
+        was added.
+
         Parameters
         ==========
         envelope_file: str
@@ -1163,19 +1171,37 @@ def from_envelope_file(envelope_file, minimum_frequency,
             Label for the names of the parameters, e.g., `recalib_H1_`
         boundary: None, 'reflective', 'periodic'
             The type of prior boundary to assign
+        correction_type: str
+            How the correction is defined, either to the :code:`data`
+            (default) or the :code:`template`. In general, data products
+            produced by the LVK calibration groups assume :code:`data`.
+            The default value will be removed in a future release and
+            this will need to be explicitly specified.
+
+            .. versionadded:: 1.4.0
 
         Returns
         =======
         prior: PriorDict
             Priors for the relevant parameters.
             This includes the frequencies of the nodes which are _not_ sampled.
         """
+        from .detector.calibration import _check_calibration_correction_type
+        correction_type = _check_calibration_correction_type(correction_type=correction_type)
+
         calibration_data = np.genfromtxt(envelope_file).T
         log_frequency_array = np.log(calibration_data[0])
-        amplitude_median = calibration_data[1] - 1
-        phase_median = calibration_data[2]
-        amplitude_sigma = (calibration_data[5] - calibration_data[3]) / 2
-        phase_sigma = (calibration_data[6] - calibration_data[4]) / 2
+
+        if correction_type.lower() == "data":
+            amplitude_median = 1 / calibration_data[1] - 1
+            phase_median = -calibration_data[2]
+            amplitude_sigma = abs(1 / calibration_data[3] - 1 / calibration_data[5]) / 2
+            phase_sigma = abs(calibration_data[6] - calibration_data[4]) / 2
+        else:
+            amplitude_median = calibration_data[1] - 1
+            phase_median = calibration_data[2]
+            amplitude_sigma = abs(calibration_data[5] - calibration_data[3]) / 2
+            phase_sigma = abs(calibration_data[6] - calibration_data[4]) / 2
 
         log_nodes = np.linspace(np.log(minimum_frequency),
                                 np.log(maximum_frequency), n_nodes)

test/gw/detector/calibration_test.py

@@ -1,5 +1,6 @@
 import os
 import unittest
+from parameterized import parameterized
 
 import numpy as np
 from bilby.core.prior import PriorDict
@@ -109,7 +110,8 @@ def test_generate_draws(self):
         self.assertEqual(draws.shape, (self.number_of_draws, sum(self.ifo.frequency_mask)))
         self.assertListEqual(list(self.priors.keys()), list(parameters.keys()))
 
-    def test_read_write_matches(self):
+    @parameterized.expand([("template",), ("data",), (None,)])
+    def test_read_write_matches(self, correction_type):
         draws, parameters = calibration._generate_calibration_draws(
             self.ifo, self.priors, self.number_of_draws
         )
@@ -119,12 +121,14 @@ def test_read_write_matches(self):
             frequency_array=frequencies,
             calibration_draws=draws,
             calibration_parameter_draws=parameters,
+            correction_type=correction_type,
         )
         self.assertTrue(os.path.exists(self.filename))
         loaded_draws, loaded_parameters = calibration.read_calibration_file(
             filename=self.filename,
             frequency_array=frequencies,
             number_of_response_curves=self.number_of_draws,
+            correction_type=correction_type,
         )
         self.assertLess(np.max(np.abs(loaded_draws - draws)), 1e-13)
         self.assertTrue(parameters.equals(loaded_parameters))