CI development

This MR will include many useful unit tests in the continuous integration. This long overdue upgrade will save a lot of manual effort in making plots for reviews, reduce the occurrence of bugs in the code, and make bugs easier to catch right away. Here are the tests that are now included:

• Run 3 versions of the calibration pipeline, and check that the right number of calibrated frames are produced
• Check the ASDs of GDS-CALIB_STRAIN and GDS-CALIB_STRAIN_CLEAN, and compare them to a "standard", nominally correct ASD. Plots are included in the CI job artifacts.
• Check the TDCFs, and compare them to nominally correct TDCFs. Plots included.
• Check the calibration state vector, and compare each bit to a nominally correct state vector. Plots included.
• Measure ratios of h(t) / Pcal and h(t) / Actuation injections, and compare to a nominally correct version. Time-series plots included.
• Use FrDiff to compare GDS-CALIB_STRAIN_CLEAN in frames produced by two otherwise identical pipelines with different start times. They are required to be identical.
• Measure the latency in a simulated low-latency-like configuration. Check that the mean latency is not too large, that there are not frequent or large excursions, and that the latency does not increase steadily with time.

Job artifacts from the most recent commit (at the time of writing) can be seen here, including plots and txt files with the data being plotted.

Getting all stages of the CI to work also required me to address some packaging issues. If I did things correctly, this MR Closes #7 (closed), #12 (closed), #14 (closed), #15 (closed), #26 (closed). I was not able to add @duncanmmacleod as a reviewer, but he created several of these issues.

tests/tests_pytest/plot.py

-#!/usr/bin/env python3
-# Copyright (C) 2018  Aaron Viets, Alexander Bartoletti
+# Copyright (C) 2023  Aaron Viets, Alexander Bartoletti
 #
 # This program is free software; you can redistribute it and/or modify it
 # under the terms of the GNU General Public License as published by the
@@ -15,175 +14,355 @@
 # with this program; if not, write to the Free Software Foundation, Inc.,
 # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
+import matplotlib; matplotlib.use('Agg')
 import numpy as np
+import time
+from math import erf
+import time
+from matplotlib import rc
 import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
+from tests.tests_pytest.ticks_and_grid import ticks_and_grid
 
 
-def plot_ASD(hoft_f='hoft_asd.txt', clean_f='clean_asd.txt'):
-    hoft_f = 'tests/tests_pytest/ASD_data/' + str(hoft_f)
-    clean_f = 'tests/tests_pytest/ASD_data/' + str(clean_f)
-    hoft = np.loadtxt(hoft_f)
-    clean = np.loadtxt(clean_f)
-
-    hoft = np.transpose(hoft)
-    clean = np.transpose(clean)
-
-    plt.rcParams["figure.figsize"] = [7.50, 7.50]
-    plt.rcParams["figure.autolayout"] = True
-    hx = hoft[0]
-    hy = hoft[1]
-    cx = clean[0]
-    cy = clean[1]
-
-#hoft
-    plt.subplot(2,1,1)
-    plt.gca().set_xscale('log')
-    plt.gca().set_yscale('log')
-    plt.title('Hoft ASD')
-    plt.plot(hx,hy, color='blue')
-    plt.xlabel('Frequency')
-    plt.ylabel('Strain h(f)')
-
-#clean
-    plt.subplot(2,1,2)
-    plt.gca().set_xscale('log')
-    plt.gca().set_yscale('log')
-    plt.title('Clean ASD')
-    plt.plot(cx,cy, color='red')
-    plt.xlabel('Frequency')
-    plt.ylabel('Strain h(f)')
-
-    plt.savefig('tests/tests_pytest/test_ASD_plot.png')
-
-
-def plot_pcal():
-    #freq_list = [102.1, 1083.3, 1083.7, 17.1, 283.9, 33.4, 410.2, 410.3, 56.4, 77.7]
-    freq_list = [17.1, 410.3, 1083.7]
-    f_list_a = []
-    f_list_e = []
-    for freq in freq_list:
-        f_list_a.append('tests/tests_pytest/pcal_data/H1_GDS-CALIB_STRAIN_over_Pcal_0_at_%sHz.txt' % freq)
-        f_list_e.append('tests/tests_pytest/pcal_data/H1_GDS-CALIB_STRAIN_over_Pcal_1_at_%sHz.txt' % freq)
-    
-    plt.rcParams["figure.figsize"] = [7.50, 10]
-    plt.rcParams["figure.autolayout"] = True
-    
-    np_lst_a = []
-    cnt = 0
-    for f in f_list_a:
-        name = freq_list[cnt]
-        cnt += 1
-        arr = np.loadtxt(f)
-        arr = np.transpose(arr)
-        x_time = arr[0]
-        y_phase = arr[2]                #magnitude and phase may be switched
-        y_magnitude = arr[1]
-
-        plt.subplot(1, 2, 1)
-        plt.title('Magnitude at Frequency %s' % name)
-        plt.plot(x_time, y_magnitude, color='green')
-        plt.xlabel('Time')
-        plt.ylabel('Magnitude')
-
-        plt.subplot(1, 2, 2)
-        plt.title('Phase at Frequency %s' % name)
-        plt.plot(x_time, y_phase, color='green')
-        plt.xlabel('Time')
-        plt.ylabel('Phase')
-        
-        plt.savefig('tests/tests_pytest/test_approx_pcal_plot_%s.png' % name)
-        plt.close()
-
-    np_lst_e = []
-    cnt = 0
-    for f in f_list_e:
-        name = freq_list[cnt]
-        cnt += 1
-        arr = np.loadtxt(f)
-        arr = np.transpose(arr)
-        x_time = arr[0]
-        y_phase = arr[2]                #magnitude and phase may be switched
-        y_magnitude = arr[1]
-
-        plt.subplot(1, 2, 1)
-        plt.title('Magnitude at Frequency %s' % name)
-        plt.plot(x_time, y_magnitude, color='green')
-        plt.xlabel('Time')
-        plt.ylabel('Magnitude')
-
-        plt.subplot(1, 2, 2)
-        plt.title('Phase at Frequency %s' % name)
-        plt.plot(x_time, y_phase, color='green')
-        plt.xlabel('Time')
-        plt.ylabel('Phase')
-
-        plt.savefig('tests/tests_pytest/test_exact_pcal_plot_%s.png' % name)
-        plt.close()
-
-
-def plot_act():
-    freq_list = [15.6, 16.4, 17.6]
-    f_list_a = []
-    f_list_e = []
-    for freq in freq_list:
-        f_list_a.append('tests/tests_pytest/act_data/H1_GDS-CALIB_STRAIN_over_Act_0_at_%sHz.txt' % freq)
-        f_list_e.append('tests/tests_pytest/act_data/H1_GDS-CALIB_STRAIN_over_Act_1_at_%sHz.txt' % freq)
-
-    plt.rcParams["figure.figsize"] = [7.50, 10]
-    plt.rcParams["figure.autolayout"] = True
-
-    np_lst_a = []
-    acnt = 0
-    for f in f_list_a:
-        name = freq_list[acnt]
-        acnt += 1
-        arr = np.loadtxt(f)
-        arr = np.transpose(arr)
-        x_time = arr[0]
-        y_phase = arr[2]                #magnitude and phase may be switched
-        y_magnitude = arr[1]
-
-        plt.subplot(1, 2, 1)
-        plt.title('Magnitude at Frequency %s' % name)
-        plt.plot(x_time, y_magnitude, color='coral')
-        plt.xlabel('Time')
-        plt.ylabel('Magnitude')
-
-
-        plt.subplot(1, 2, 2)
-        plt.title('Phase at Frequency %s' % name)
-        plt.plot(x_time, y_phase, color='coral')
-        plt.xlabel('Time')
-        plt.ylabel('Phase')
-
-        plt.savefig('tests/tests_pytest/test_approx_act_plot_%s.png' % name)
-        plt.close()
-
-    np_lst_e = []
-    ecnt = 0
-    for f in f_list_e:
-        name = freq_list[ecnt]
-        ecnt += 1
-        arr = np.loadtxt(f)
-        arr = np.transpose(arr)
-        x_time = arr[0]
-        y_phase = arr[2]                #magnitude and phase may be switched
-        y_magnitude = arr[1]
-
-        plt.subplot(1, 2, 1)
-        plt.title('Magnitude at Frequency %s' % name)
-        plt.plot(x_time, y_magnitude, color='coral')
-        plt.xlabel('Time')
-        plt.ylabel('Magnitude')
-
-
-        plt.subplot(1, 2, 2)
-        plt.title('Phase at Frequency %s' % name)
-        plt.plot(x_time, y_phase, color='coral')
-        plt.xlabel('Time')
-        plt.ylabel('Phase')
-
-        plt.savefig('tests/tests_pytest/test_exact_act_plot_%s.png' % name)
-        plt.close()
+# Function to compute standard deviation that is not strongly affected by outliers
+def stdmed(x):
+	x = np.asarray(x)
+	N = len(x)
+	xmed = np.median(x)
+	diff = x - xmed
+	square_diff = diff * diff
+	diff_med = np.percentile(square_diff, 100 * erf(1 / np.sqrt(2)))
+	return np.sqrt(diff_med * N / (N - 1))
 
 
+def plot_ASD(hoft_f = 'Approx_hoft_asd.txt', clean_f = 'Approx_clean_asd.txt', standard = None):
+	version = hoft_f.split('_')[0]
+	plot_filename = 'tests/tests_pytest/ASD_data/' + hoft_f.split('.')[0] + '_plot.png'
+	hoft_f = 'tests/tests_pytest/ASD_data/' + str(hoft_f)
+	clean_f = 'tests/tests_pytest/ASD_data/' + str(clean_f)
+	labels = ['strain', 'clean']
+	colors = ['red', 'limegreen']
+	hoft = np.loadtxt(hoft_f)
+	clean = np.loadtxt(clean_f)
+
+	hoft = np.transpose(hoft)
+	clean = np.transpose(clean)
+
+	matplotlib.rcParams['font.size'] = 22
+	matplotlib.rcParams['legend.fontsize'] = 18
+	plt.rcParams["figure.figsize"] = [10, 6]
+	# plt.rcParams["figure.autolayout"] = True
+	hx = hoft[0]
+	hy = hoft[1]
+	cx = clean[0]
+	cy = clean[1]
+
+	if standard is not None:
+		standard = np.transpose(np.loadtxt('tests/tests_pytest/ASD_data/' + standard))[1]
+		hy /= standard
+		cy /= standard
+		plot_filename = plot_filename.split('.')[0] + "_ratio.png"
+	# hoft
+	plt.title('Strain ASD with %s TDCFs' % version)
+	plt.plot(hx, hy, color = colors[0], linewidth = 0.75)
+	patches = [mpatches.Patch(color = colors[j], label = labels[j]) for j in range(len(labels))]
+	plt.legend(handles = patches, loc = 'upper right', ncol = 1)
+	if standard is None:
+		plt.title('Strain ASD with %s TDCFs' % version)
+		plt.ylabel(r'${\rm ASD}\ \left[{\rm strain / }\sqrt{\rm Hz}\right]$')
+		ticks_and_grid(plt.gca(), xmin = 1, xmax = 8192, ymin = 1e-24, ymax = 1e-17, xscale = 'log', yscale = 'log')
+	else:
+		plt.title('Strain ASD / Standard with %s TDCFs' % version)
+		plt.ylabel('ASD Ratio')
+		ticks_and_grid(plt.gca(), xmin = 1, xmax = 8192, ymin = 0.8, ymax = 1.2, xscale = 'log', yscale = 'linear')
+	plt.xlabel(r'${\rm Frequency \ [Hz]}$')
+	plt.tight_layout()
+
+	# clean
+	plt.plot(cx,cy, color = colors[1], linewidth = 0.75)
+
+	plt.savefig(plot_filename)
+	plt.close()
+
+def plot_TDCFs():
+	matplotlib.rcParams['font.size'] = 20
+	matplotlib.rcParams['legend.fontsize'] = 20
+
+	ifo = 'H1'
+	labels = ['Approx', 'Exact']
+	colors = ['red', 'limegreen']
+	TDCF_names = ['KAPPA_TST_REAL', 'KAPPA_TST_IMAGINARY', 'KAPPA_PUM_REAL', 'KAPPA_PUM_IMAGINARY', 'KAPPA_UIM_REAL', 'KAPPA_UIM_IMAGINARY', 'KAPPA_C', 'F_CC', 'F_S_SQUARED', 'SRC_Q_INVERSE']
+	plot_min = [0.9, -0.1, 0.9, -0.1, 0.9, -0.1, 0.9, 400, -100, -2]
+	plot_max = [1.1, 0.1, 1.1, 0.1, 1.1, 0.1, 1.1, 500, 100, 2]
+
+	for i in range(len(TDCF_names)):
+		Approx_data = np.transpose(np.loadtxt("tests/tests_pytest/TDCF_data/%s_Approx.txt" % TDCF_names[i]))
+		Exact_data = np.transpose(np.loadtxt("tests/tests_pytest/TDCF_data/%s_Exact.txt" % TDCF_names[i]))
+		Approx_t = Approx_data[0]
+		Exact_t = Exact_data[0]
+		Approx_TDCF = Approx_data[1]
+		Exact_TDCF = Exact_data[1]
+
+		t_start = Approx_t[0]
+		dur = Approx_t[-1] - t_start
+		t_unit = 'seconds'
+		sec_per_t_unit = 1.0
+		if dur > 60 * 60 * 100:
+			t_unit = 'days'
+			sec_per_t_unit = 60.0 * 60.0 * 24.0
+		elif dur > 60 * 100:
+			t_unit = 'hours'
+			sec_per_t_unit = 60.0 * 60.0
+		elif dur > 100:
+			t_unit = 'minutes'
+			sec_per_t_unit = 60.0
+
+		Approx_t -= t_start
+		Exact_t -= t_start
+		Approx_t /= sec_per_t_unit
+		Exact_t /= sec_per_t_unit
+
+		# Make plots
+		plt.figure(figsize = (10, 6))
+		plt.plot(Approx_t, Approx_TDCF, color = 'red', linewidth = 0.75, label = 'Approx')
+		plt.plot(Exact_t, Exact_TDCF, color = 'limegreen', linewidth = 0.75, label = 'Exact')
+		patches = [mpatches.Patch(color = colors[j], label = labels[j]) for j in range(len(labels))]
+		plt.legend(handles = patches, loc = 'upper right', ncol = 1)
+		plt.xlabel("Time in %s since %s UTC" % (t_unit, time.strftime("%b %d %Y %H:%M:%S", time.gmtime(t_start + 315964782))))
+		if "F_S_SQUARED" in TDCF_names[i]:
+			plt.ylabel(r"%s [${\rm Hz}^2$]" % TDCF_names[i])
+		elif "F_CC" in TDCF_names[i]:
+			plt.ylabel("%s [Hz]" % TDCF_names[i])
+		else:
+			plt.ylabel(TDCF_names[i])
+
+		ticks_and_grid(plt.gca(), ymin = plot_min[i], ymax = plot_max[i])
+		plt.tight_layout()
+
+		plt.savefig("tests/tests_pytest/TDCF_data/%s.png" % TDCF_names[i])
+		plt.close()
+
+def plot_deltal_over_inj_timeseries(inj):
+	ifo = 'H1'
+	plot_labels = ['Approx','Exact']
+	channels = ['GDS-CALIB_STRAIN', 'GDS-CALIB_STRAIN']
+	coherence_time = 148
+	if inj == 'act':
+		frequencies = [17.6, 16.4, 15.6] # Hz
+		inj_names = ['TST_exc', 'PUM_exc', 'UIM_exc']
+		stages = ['T', 'P', 'U']
+	elif inj == 'pcaly':
+		frequencies = [17.1, 410.3, 1083.7] # Hz
+		inj_names = ['Pcal', 'Pcal', 'Pcal']
+		stages = ['pcy', 'pcy', 'pcy']
+	elif inj == 'pcalx':
+		frequencies = [33.43, 77.73, 102.13, 283.91] # Hz
+		inj_names = ['Pcal', 'Pcal', 'Pcal', 'Pcal']
+		stages = ['pcx', 'pcx', 'pcx', 'pcx']
+	else:
+		assert False
+
+	matplotlib.rcParams['font.size'] = 32
+	matplotlib.rcParams['legend.fontsize'] = 16
+
+	# Read data from files and plot it
+	colors = ['red', 'limegreen', 'mediumblue', 'gold', 'b', 'm'] # Hopefully the user will not want to plot more than six datasets on one plot.
+	for i in range(0, len(frequencies)):
+		for j in range(len(channels)):
+			data = np.loadtxt("tests/tests_pytest/%s_data/%s_%s_over_%s_%d_at_%0.1fHz.txt" % (inj[:4], ifo, channels[j], inj_names[i], j, frequencies[i]))
+			if j == 0:
+				t_start = data[0][0]
+				dur = data[len(data) - 1][0] - t_start
+				dur_in_seconds = dur
+				t_unit = 'seconds'
+				sec_per_t_unit = 1.0
+				if dur > 60 * 60 * 100:
+					t_unit = 'days'
+					sec_per_t_unit = 60.0 * 60.0 * 24.0
+				elif dur > 60 * 100:
+					t_unit = 'hours'
+					sec_per_t_unit = 60.0 * 60.0
+				elif dur > 100:
+					t_unit = 'minutes'
+					sec_per_t_unit = 60.0
+
+				markersize = 150.0 * np.sqrt(float(coherence_time / dur))
+				markersize = min(markersize, 10.0)
+				markersize = max(markersize, 0.2)
+
+				times = []
+				magnitudes = []
+				phases = []
+
+			times.append([])
+			magnitudes.append([])
+			phases.append([])
+			for k in range(len(data)):
+				if(data[k][0] % coherence_time == 0):
+					times[j].append((data[k][0] - t_start) / sec_per_t_unit)
+					magnitudes[j].append(data[k][1])
+					phases[j].append(data[k][2])
+			# Make plots
+			if i == 0 and j == 0:
+				plt.figure(figsize = (8 * len(frequencies), 15))
+			plt.subplot(2, len(frequencies), i + 1)
+			plt.plot(times[j], magnitudes[j], colors[j], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s} \ [\mu_{1/2} = %0.4f, \sigma = %0.4f]$' % (plot_labels[j], np.median(magnitudes[j]), stdmed(magnitudes[j])))
+			if j == 0:
+				plt.title(r'${\rm %s} \ \widetilde{\Delta L}_{\rm free} / \tilde{x}_{\rm %s} \  {\rm at \  %0.1f \  Hz}$' % (ifo, stages[i], frequencies[i]), fontsize = 32)
+				ticks_and_grid(plt.gca(), ymin = 0.9, ymax = 1.1)
+				if i == 0:
+					plt.ylabel(r'${\rm Magnitude}$')
+			leg = plt.legend(fancybox = True, markerscale = 16.0 / markersize, numpoints = 1, loc = 'upper right')
+			leg.get_frame().set_alpha(0.8)
+			plt.subplot(2, len(frequencies), len(frequencies) + i + 1)
+			plt.plot(times[j], phases[j], colors[j], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s} \ [\mu_{1/2} = %0.3f^{\circ}, \sigma = %0.3f^{\circ}]$' % (plot_labels[j], np.median(phases[j]), stdmed(phases[j])))
+			leg = plt.legend(fancybox = True, markerscale = 16.0 / markersize, numpoints = 1, loc = 'upper right')
+			leg.get_frame().set_alpha(0.8)
+			if j == 0:
+				ticks_and_grid(plt.gca(), ymin = -6, ymax = 6)
+				if len(frequencies) < 3 or i == int((len(frequencies) - 0.1) / 2.0):
+					plt.xlabel("Time in %s since %s UTC" % (t_unit, time.strftime("%b %d %Y %H:%M:%S", time.gmtime(t_start + 315964782))))
+				if i == 0:
+					plt.ylabel(r'${\rm Phase \  [deg]}$')
+	plt.savefig("tests/tests_pytest/%s_data/%s_deltal_over_%s_%d-%d.png" % (inj[:4], ifo, inj, int(t_start), int(dur_in_seconds)))
+	plt.close()
+
+def plot_statevector(filename):
+	matplotlib.rcParams['font.size'] = 20
+	matplotlib.rcParams['legend.fontsize'] = 20
+
+	ifo = 'H1'
+
+	# Bit definitions
+	bits = []
+	bits.append("hoft_ok")
+	bits.append("obs_intent")
+	bits.append("lownoise")
+	bits.append("filters_ok")
+	bits.append("no_gap")
+	bits.append("no_stoch_inj")
+	bits.append("no_cbc_inj")
+	bits.append("no_burst_inj")
+	bits.append("no_detchar_inj")
+	bits.append("undisturbed_ok")
+	bits.append("ktst_smooth")
+	bits.append("kpum_smooth")
+	bits.append("kuim_smooth")
+	bits.append("kc_smooth")
+	bits.append("fcc_smooth")
+	bits.append("fs_smooth")
+	bits.append("fs_over_Q_smooth")
+	bits.append("line_sub")
+	bits.append("noise_sub")
+	bits.append("noise_sub_gate")
+	bits.append("nonsens_sub")
+
+	N_bits = len(bits)
+
+	statevector = np.loadtxt('tests/tests_pytest/State_Vector_data/%s.txt' % filename)
+
+	t = np.transpose(statevector)[0]
+	t_start = t[0]
+	t -= t_start
+	sr = 16
+	cadence = 1 # Number of state vector samples to combine into one sample on the plot
+	while len(t) / cadence / 2 > 500:
+		cadence *= 2
+	width = cadence / sr
+	t = t[::cadence]
+	dur = t[-1]
+	t_unit = 'seconds'
+	sec_per_t_unit = 1.0
+	if dur > 60 * 60 * 100:
+		t_unit = 'days'
+		sec_per_t_unit = 60.0 * 60.0 * 24.0
+	elif dur > 60 * 100:
+		t_unit = 'hours'
+		sec_per_t_unit = 60.0 * 60.0
+	elif dur > 100:
+		t_unit = 'minutes'
+		sec_per_t_unit = 60.0
+	t /= sec_per_t_unit
+
+	t5 = int(round(t[-1] / 5))
+	tier = int(np.round(3*np.log10(t5)))
+	t_tick_space = 10**(tier // 3) * int(round(10**((tier % 3) / 3)))
+	t_ticks = np.arange(0, t[-1], t_tick_space)
+	t_tick_locations = t_ticks * sec_per_t_unit
+
+	statevector = np.transpose(statevector)[1]
+	bit_values = np.ones((N_bits, len(t)), dtype = int)
+	for i in range(len(statevector)):
+		sv = bin(int(round(statevector[i])))
+		if len(sv) - 2 > N_bits:
+			raise ValueError("Found %d bits, but only expected %d bits" % (len(sv) - 2, N_bits))
+		for j in range(len(sv) - 2):
+			bit_values[j][i // cadence] &= int(sv[len(sv) - 1 - j])
+		# The rest of the bits are zero
+		for j in range(len(sv) - 2, N_bits):
+			bit_values[j][i // cadence] = 0
+
+	# Color of the bars on the state vector plot
+	on_color = 'green'
+	off_color = 'red'
+
+	# Make plots
+	fig, axs = plt.subplots(N_bits, figsize = (15, 8))
+	fig.suptitle("%s Calibration State" % ifo)
+
+	for i in range(N_bits):
+		colors = [(on_color if j > 0.9 else off_color) for j in bit_values[i]]
+
+		for j, k in enumerate(bit_values[i]):
+			axs[i].barh(y = 0, width = width, height = (k + 2) / 3, left = j * width, color = colors[j])
+		axs[i].barh(y = 0, width = 1, height = 1, left = len(bit_values[i]) * width, color = 'black')
+		plt.sca(axs[i])
+		if i < N_bits - 1:
+			plt.xticks(ticks = t_tick_locations)
+			axs[i].xaxis.set_tick_params(labelbottom=False)
+		axs[i].yaxis.set_tick_params(left=False)
+		plt.yticks(ticks = [0, 0.25, 0.5], labels = [bits[i], '', ''])
+
+	plt.xticks(ticks = t_tick_locations, labels = t_ticks)
+	plt.xlabel("Time in %s since %s UTC" % (t_unit, time.strftime("%b %d %Y %H:%M:%S", time.gmtime(t_start + 315964782))))
+
+	#fig.tight_layout(rect = [0.2, 0, 1, 1])
+	#fig.tight_layout()
+	fig.subplots_adjust(left = 0.2)
+
+	plt.savefig('tests/tests_pytest/State_Vector_data/%s_%s_plot_%d-%d.png' % (ifo, filename, t_start, dur))
+	plt.close()
+
+def plot_latency(o_gps, lat):
+	t0 = o_gps[0]
+	t = o_gps - t0
+	dur = t[-1]
+	dur_in_seconds = dur
+	t_unit = 'seconds'
+	sec_per_t_unit = 1.0
+	if dur > 60 * 60 * 100:
+		t_unit = 'days'
+		sec_per_t_unit = 60.0 * 60.0 * 24.0
+	elif dur > 60 * 100:
+		t_unit = 'hours'
+		sec_per_t_unit = 60.0 * 60.0
+	elif dur > 100:
+		t_unit = 'minutes'
+		sec_per_t_unit = 60.0
+	t /= sec_per_t_unit
+
+	matplotlib.rcParams['font.size'] = 22
+	matplotlib.rcParams['legend.fontsize'] = 18
+	plt.rcParams["figure.figsize"] = [12, 6]
+
+	plt.plot(t, lat, linewidth = 0.75)
+	plt.xlabel("Time in %s since %s UTC" % (t_unit, time.strftime("%b %d %Y %H:%M:%S", time.gmtime(int(t0) + 315964782))))
+	plt.ylabel("Latency [s]")
+	plt.title("gstlal calibration pipeline latency")
+	ticks_and_grid(plt.gca(), ymin = -1, ymax = 10)
+	plt.tight_layout()
+
+	plt.savefig("tests/tests_pytest/gstlal_cal_latency.png")
+	plt.close()
+