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

@@ -231,9 +231,9 @@ class CubicSpline(Recalibrate):
         log10f_per_deltalog10f = (
             np.log10(frequency_array) - self.log_spline_points[0]
         ) / self._delta_log_spline_points
-        previous_idxs = np.floor(log10f_per_deltalog10f).astype(int)
-        next_idxs = previous_idxs + 1
-        b = log10f_per_deltalog10f - previous_idxs
+        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
@@ -241,13 +241,13 @@ class CubicSpline(Recalibrate):
         self.set_calibration_parameters(**params)
         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_idxs] + b * amplitude_parameters[next_idxs] + \
-            c * _spline_coefficients[previous_idxs] + d * _spline_coefficients[next_idxs]
+        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_idxs] + b * phase_parameters[next_idxs] + \
-            c * _spline_coefficients[previous_idxs] + d * _spline_coefficients[next_idxs]
+        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)