Overall response change plot as first plot in Reports

I've put together a plot of how the response function of the current reports' measurements differs from the exported model: both as a raw difference, as well as the kappa corrected difference.

The kappa corrected difference can be used to quickly see if there is significant systematic error between your current measurements and the running calibibration pipeline after kappa corrections.

Attachment is what it looks like:

pydarm/cmd/_history.py

 from copy import deepcopy
 
 import numpy as np
+from gwpy.timeseries import TimeSeriesDict as tsd
 from scipy.signal import freqresp
 from matplotlib import pyplot as plt
 from matplotlib.ticker import MultipleLocator
+from joblib import Memory
 
 from ._log import CMDError
 from ._const import FREE_PARAM_LABEL_MAP
@@ -42,6 +44,130 @@ def adjust_mag_yticks(ax, maxsteps=7):
     ax.yaxis.set_minor_locator(loc_minor)
 
 
+def response_history(new_report, previous_exported_report):
+    """generate plot of overall response chane
+    in 2 flavours: raw response change and kappa corrected response change
+    """
+
+    # fetch 1 hour average of kappas:
+    @Memory('_cache').cache
+    def get_kappas(IFO, gps_time):
+        '''fetch kappas from NDS server
+            implicit assumptions being made here:
+               -IFO is well thermalised when you make the measurements
+               -you have at least 1 hours' worth of kappa data which you can draw a good average
+            1 hour comes from testing in LLO alog 72914, to make measurements made in identical
+               configuration reliably give the same response change
+        '''
+        end = int(gps_time)-300  # go back 5 minutes in case you run broadband before sweeps
+        start = end - 3600
+        kappa_data = tsd.get(
+            [f'{IFO}:CAL-CS_TDEP_KAPPA_C_OUTPUT',
+             f'{IFO}:CAL-CS_TDEP_F_C_OUTPUT',
+             f'{IFO}:CAL-CS_TDEP_KAPPA_UIM_OUTPUT',
+             f'{IFO}:CAL-CS_TDEP_KAPPA_PUM_OUTPUT',
+             f'{IFO}:CAL-CS_TDEP_KAPPA_TST_OUTPUT'],
+            start, end,
+            allow_tape=True,
+            verbose=False,
+        )
+
+        kappa_dict = {'c': np.mean(kappa_data[f'{IFO}:CAL-CS_TDEP_KAPPA_C_OUTPUT'].data),
+                      'f_c': np.mean(kappa_data[f'{IFO}:CAL-CS_TDEP_F_C_OUTPUT'].data),
+                      'uim': np.mean(kappa_data[f'{IFO}:CAL-CS_TDEP_KAPPA_UIM_OUTPUT'].data),
+                      'pum': np.mean(kappa_data[f'{IFO}:CAL-CS_TDEP_KAPPA_PUM_OUTPUT'].data),
+                      'tst': np.mean(kappa_data[f'{IFO}:CAL-CS_TDEP_KAPPA_TST_OUTPUT'].data)}
+        return kappa_dict
+
+    def apply_kappa_to_model(kappas, model):
+        model.sensing.coupled_cavity_optical_gain *= kappas['c']
+        model.sensing.coupled_cavity_pole_frequency = kappas['f_c']
+        # any of the below can fail if the actuation stage is not being used in the ifo/ini
+        try:
+            model.actuation.xarm.uim_npa *= kappas['uim']
+        except Exception:
+            pass
+        try:
+            model.actuation.yarm.uim_npa *= kappas['uim']
+        except Exception:
+            pass
+        try:
+            model.actuation.xarm.pum_npa *= kappas['pum']
+        except Exception:
+            pass
+        try:
+            model.actuation.yarm.pum_npa *= kappas['pum']
+        except Exception:
+            pass
+        try:
+            model.actuation.xarm.tst_npv2 *= kappas['tst']
+        except Exception:
+            pass
+        try:
+            model.actuation.yarm.tst_npv2 *= kappas['tst']
+        except Exception:
+            pass
+        return model
+
+    new_report = deepcopy(new_report)
+    old_report = deepcopy(previous_exported_report)
+    # apply kappas to the old model to propagate the state of the currently running (exported)
+    # model relative to new measurements
+    old_report_kappa_correct_model = apply_kappa_to_model(
+                                        get_kappas(new_report.IFO, new_report.gps()),
+                                        old_report.model)
+
+    # prepare frequency responses
+    freq = np.geomspace(1, 8000, 100000)
+    R_new = new_report.model.compute_response_function(freq)
+    R_old = old_report.model.compute_response_function(freq)
+    R_old_kappa_correct = old_report_kappa_correct_model.compute_response_function(freq)
+
+    # main figure
+    plt.rcParams.update({'mathtext.default': 'regular'})
+    resp_hist_fig, axs = plt.subplots(2, 1)
+    resp_hist_fig.suptitle(
+        'Response change since the last exported model\n',
+        fontsize=20)
+    s = f'Responses drawn from {new_report.model_file} relative to export {old_report.model_file}'
+    resp_hist_fig.text(
+        .5, .93,
+        s,
+        horizontalalignment='center',
+        transform=resp_hist_fig.transFigure)
+    resp_hist_fig.tight_layout(rect=(0.05, 0.05, 1, 1))
+
+    # plot parameters
+    freq_plot_lims = [7.5, 5000]
+    mag_ratio_plot_lims = [0.95, 1.05]
+    phase_ratio_plot_lims = [-5, 5]
+
+    axs[0].semilogx(freq, np.abs(R_new/R_old), color='blue', linestyle='--', linewidth=0.5,
+                    label='R_current_model / R_last_exported_model [Raw]')
+    axs[0].semilogx(freq, np.abs(R_new/R_old_kappa_correct), color='red',
+                    label='R_current_model / R_last_exported_model [Kappa Corrected]')
+    axs[0].legend()
+    axs[0].set(ylabel=r'$R_{pyDARM} / R_{pyDARM_0}$ [mag]')
+    axs[0].set_xlim(freq_plot_lims[0], freq_plot_lims[1])
+    axs[0].set_ylim(mag_ratio_plot_lims[0], mag_ratio_plot_lims[1])
+    axs[0].set_yticks(np.arange(mag_ratio_plot_lims[0], mag_ratio_plot_lims[1], 0.01))
+
+    axs[1].semilogx(freq, np.angle(R_new/R_old, deg=True), color='b', linestyle='--', linewidth=0.5)
+    axs[1].semilogx(freq, np.angle(R_new/R_old_kappa_correct, deg=True), color='red')
+    axs[1].set(xlabel='Freq (Hz)')
+    axs[1].set(ylabel=r'$R_{pyDARM} / R_{pyDARM_0}$ [deg]')
+    axs[1].set_xlim(7.5, 5000)
+    axs[1].set_ylim(phase_ratio_plot_lims[0], phase_ratio_plot_lims[1])
+    axs[1].set_yticks(np.arange(phase_ratio_plot_lims[0], phase_ratio_plot_lims[1]+1, 1))
+    # Wrap up and save figure
+    resp_hist_fig.savefig(
+        new_report.gen_path(
+            'response_history.png'
+        )
+    )
+    return resp_hist_fig
+
+
 def sensing_history(new_report, epoch_reports, ref_report):
     """generate sensing history plots