REFACTOR: make optimized spline calibration backward compatible

bilby/gw/detector/calibration.py

@@ -181,24 +181,38 @@ 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))
+    @property
+    def delta_log_spline_points(self):
+        if not hasattr(self, "__delta_log_spline_points"):
+            self._delta_log_spline_points = self._log_spline_points[1] - self._log_spline_points[0]
+        return self._delta_log_spline_points
+
+    @property
+    def nodes_to_spline_coefficients(self):
+        if not hasattr(self, "_nodes_to_spline_coefficients"):
+            self._setup_spline_coefficients()
+        return self._nodes_to_spline_coefficients
+
+    def _setup_spline_coefficients(self):
+        """
+        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=(self.n_points, self.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):
+        for i in range(1, self.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))
+        tmp2 = np.zeros(shape=(self.n_points, self.n_points))
         for i in range(1, n_points - 1):
             tmp2[i, i - 1] = 1
             tmp2[i, i] = -2
@@ -213,6 +227,18 @@ class CubicSpline(Recalibrate):
         return self.__class__.__name__ + '(prefix=\'{}\', minimum_frequency={}, maximum_frequency={}, n_points={})'\
             .format(self.prefix, self.minimum_frequency, self.maximum_frequency, self.n_points)
 
+    def _evaluate_spline(self, kind, a, b, c, d, previous_nodes):
+        """Evaluate Eq. (1) in https://dcc.ligo.org/LIGO-T2300140"""
+        parameters = np.array([self.params[f"{kind}_{ii}"] for ii in range(self.n_points)])
+        next_nodes = previous_nodes + 1
+        spline_coefficients = self.nodes_to_spline_coefficients.dot(amplitude_parameters)
+        return (
+            a * amplitude_parameters[previous_nodes]
+            + b * amplitude_parameters[next_nodes]
+            + c * spline_coefficients[previous_nodes]
+            + d * spline_coefficients[next_nodes]
+        )
+
     def get_calibration_factor(self, frequency_array, **params):
         """Apply calibration model
 
@@ -233,25 +259,17 @@ class CubicSpline(Recalibrate):
         """
         log10f_per_deltalog10f = (
             np.log10(frequency_array) - self.log_spline_points[0]
-        ) / self._delta_log_spline_points
+        ) / 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 = 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]
 
+        delta_amplitude = self._evaluate_spline("amplitude", a, b, c, d, previous_nodes)
+        delta_phase = self._evaluate_spline("phase", a, b, c, d, previous_nodes)
         calibration_factor = (1 + delta_amplitude) * (2 + 1j * delta_phase) / (2 - 1j * delta_phase)
 
         return calibration_factor