bilby/gw/detector/calibration.py: extrapolate calibration errors for smooth transition

This fixes #637 (closed). Since this is a corner case, I am happy to close this if extrapolation is not preferred for other reasons. This also benefits the multiband likelihood, as it requests waveform values beyond the high-frequency cutoff and assumes smoothness of waveform there, although log-likelihood errors introduced by the abrupt termination of calibration errors are typically below unity.

bilby/gw/detector/calibration.py

@@ -53,11 +53,11 @@ def read_calibration_file(filename, frequency_array, number_of_response_curves,
         calibration_phase = calibration_phase.view(np.float64).reshape(calibration_phase.shape + (-1,))
         calibration_frequencies = calibration_frequencies.view(np.float64)
 
-    # interpolate to the frequency array (where if outside the range of the calibration uncertainty its fixed to 1)
+    # interpolate to the frequency array (where if outside the range of the calibration uncertainty it is extrapolated)
     calibration_draws = calibration_amplitude * np.exp(1j * calibration_phase)
     calibration_draws = interp1d(
         calibration_frequencies, calibration_draws, kind='cubic',
-        bounds_error=False, fill_value=1)(frequency_array)
+        bounds_error=False, fill_value='extrapolate')(frequency_array)
 
     try:
         parameter_draws = pd.read_hdf(filename, key="CalParams")
@@ -213,13 +213,13 @@ class CubicSpline(Recalibrate):
                                 for ii in range(self.n_points)]
         delta_amplitude = interp1d(
             self.log_spline_points, amplitude_parameters, kind='cubic',
-            bounds_error=False, fill_value=0)(np.log10(frequency_array))
+            bounds_error=False, fill_value='extrapolate')(np.log10(frequency_array))
 
         phase_parameters = [
             self.params['phase_{}'.format(ii)] for ii in range(self.n_points)]
         delta_phase = interp1d(
             self.log_spline_points, phase_parameters, kind='cubic',
-            bounds_error=False, fill_value=0)(np.log10(frequency_array))
+            bounds_error=False, fill_value='extrapolate')(np.log10(frequency_array))
 
         calibration_factor = (1 + delta_amplitude) * (2 + 1j * delta_phase) / (2 - 1j * delta_phase)