From 0474f80e86302b6f955f0d64fbe7bac6d8780bb2 Mon Sep 17 00:00:00 2001 From: Soichiro Morisaki <soichiro.morisaki@ligo.org> Date: Fri, 24 Mar 2023 23:16:38 +0900 Subject: [PATCH] bilby/gw/detector/calibration.py: optimize spline interpolation of calibration uncertainties --- bilby/gw/detector/calibration.py | 53 +++++++++++++++++++++++++------- 1 file changed, 42 insertions(+), 11 deletions(-) diff --git a/bilby/gw/detector/calibration.py b/bilby/gw/detector/calibration.py index a5811c1ae..4f9aaf68d 100644 --- a/bilby/gw/detector/calibration.py +++ b/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) -- GitLab