bilby/gw/detector/calibration.py: optimize spline interpolation of calibration uncertainties

bilby/gw/detector/calibration.py

@@ -181,6 +181,29 @@ class CubicSpline(Recalibrate):
         self.maximum_frequency = maximum_frequency
         self._log_spline_points = np.linspace(
             np.log10(minimum_frequency), np.log10(maximum_frequency), n_points)
+        self._delta_log_spline_points = self._log_spline_points[1] - self._log_spline_points[0]
+
+        # Precompute matrix converting values at nodes to spline coefficients.
+        # The algorithm for interpolation is described in
+        # https://dcc.ligo.org/LIGO-T2300140, and the matrix calculated here is
+        # to solve Eq. (9) in the note.
+        tmp1 = np.zeros(shape=(n_points, n_points))
+        tmp1[0, 0] = -1
+        tmp1[0, 1] = 2
+        tmp1[0, 2] = -1
+        tmp1[-1, -3] = -1
+        tmp1[-1, -2] = 2
+        tmp1[-1, -1] = -1
+        for i in range(1, n_points - 1):
+            tmp1[i, i - 1] = 1 / 6
+            tmp1[i, i] = 2 / 3
+            tmp1[i, i + 1] = 1 / 6
+        tmp2 = np.zeros(shape=(n_points, n_points))
+        for i in range(1, n_points - 1):
+            tmp2[i, i - 1] = 1
+            tmp2[i, i] = -2
+            tmp2[i, i + 1] = 1
+        self._nodes_to_spline_coefficients = np.linalg.solve(tmp1, tmp2)
 
     @property
     def log_spline_points(self):
@@ -208,18 +231,26 @@ class CubicSpline(Recalibrate):
         calibration_factor : array-like
             The factor to multiply the strain by.
         """
+        log10f_per_deltalog10f = (
+            np.log10(frequency_array) - self.log_spline_points[0]
+        ) / self._delta_log_spline_points
+        previous_nodes = np.clip(np.floor(log10f_per_deltalog10f).astype(int), a_min=0, a_max=self.n_points - 2)
+        next_nodes = previous_nodes + 1
+        b = log10f_per_deltalog10f - previous_nodes
+        a = 1 - b
+        c = (a**3 - a) / 6
+        d = (b**3 - b) / 6
+
         self.set_calibration_parameters(**params)
-        amplitude_parameters = [self.params['amplitude_{}'.format(ii)]
-                                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))
-
-        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))
+        amplitude_parameters = np.array([self.params['amplitude_{}'.format(ii)] for ii in range(self.n_points)])
+        _spline_coefficients = self._nodes_to_spline_coefficients.dot(amplitude_parameters)
+        delta_amplitude = a * amplitude_parameters[previous_nodes] + b * amplitude_parameters[next_nodes] + \
+            c * _spline_coefficients[previous_nodes] + d * _spline_coefficients[next_nodes]
+
+        phase_parameters = np.array([self.params['phase_{}'.format(ii)] for ii in range(self.n_points)])
+        _spline_coefficients = self._nodes_to_spline_coefficients.dot(phase_parameters)
+        delta_phase = a * phase_parameters[previous_nodes] + b * phase_parameters[next_nodes] + \
+            c * _spline_coefficients[previous_nodes] + d * _spline_coefficients[next_nodes]
 
         calibration_factor = (1 + delta_amplitude) * (2 + 1j * delta_phase) / (2 - 1j * delta_phase)