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

Optimize spline interpolation of calibration uncertainties utilizing formula applicable when spline nodes are evenly spaced. This is similar to the improvement made for ROQ weights in !971 (merged).

Below is a test script used for testing this change:

from bilby.gw.detector.calibration import CubicSpline
import numpy as np

n_points = 10
flow, fhigh = 20., 2048.
prefix = "recalib_H1_"

frequencies = np.linspace(flow, fhigh, 100)

parameters = {}
for i in range(n_points):
    parameters[f"{prefix}amplitude_{i}"] = np.random.uniform(-0.1, 0.1)
    parameters[f"{prefix}phase_{i}"] = np.random.uniform(-0.1 * np.pi, 0.1 * np.pi)

calobj1 = CubicSpline(prefix, flow, fhigh, n_points)
calobj2 = CubicSplineNew(prefix, flow, fhigh, n_points)

Below are results:

In:
np.max(np.abs(
    calobj1.get_calibration_factor(frequencies, **parameters) -
    calobj2.get_calibration_factor(frequencies, **parameters)
))

Out:
3.619316683849706e-15  

In: 
%%timeit -n 100
calobj1.get_calibration_factor(frequencies, **parameters)  

Out: 
585 µs ± 106 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)  

In:
%%timeit -n 100
calobj2.get_calibration_factor(frequencies, **parameters)  

Out:  
101 µs ± 23.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)  

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)