diff --git a/bin/gstlal_compute_strain b/bin/gstlal_compute_strain
index b600ff0abbd91ab8f2c059af7cdfb59fefeb1e70..fe5311ea74a6581c6f5f09c4bd3e117cb694f220 100755
--- a/bin/gstlal_compute_strain
+++ b/bin/gstlal_compute_strain
@@ -1046,7 +1046,8 @@ if compute_calib_statevector:
noisesub_gate_channel = ChannelNames["noisesubgatechannel"] if "noisesubgatechannel" in ChannelNames else lownoise_channel_name
noisesub_gate_bitmask = int(Bitmasks["noisesubgatebitmask"]) if "noisesubgatebitmask" in Bitmasks else lownoise_bitmask
noisesub_gate_threshold = (int(Bitmasks["noisesubgatethreshold"]) if Bitmasks["noisesubgatethreshold"] != "None" else None) if "noisesubgatethreshold" in Bitmasks else None
-if compute_calib_statevector and (any(line_witness_channel_list) or any(witness_channel_list)) and noisesub_gate_channel != obsintent_channel_name and noisesub_gate_channel != lownoise_channel_name and noisesub_gate_channel != hwinj_channel_name and noisesub_gate_channel != filterclock_channel_name and noisesub_gate_channel != undisturbed_ok_channel and (noisesub_gate_bitmask > 0 or noisesub_gate_threshold is not None) and noisesub_gate_channel != "None" and not noisesub_use_filter_clock:
+ noisesub_gate_value = (int(Bitmasks["noisesubgatevalue"]) if Bitmasks["noisesubgatevalue"] != "None" else None) if "noisesubgatevalue" in Bitmasks else None
+if compute_calib_statevector and (any(line_witness_channel_list) or any(witness_channel_list)) and noisesub_gate_channel != obsintent_channel_name and noisesub_gate_channel != lownoise_channel_name and noisesub_gate_channel != hwinj_channel_name and noisesub_gate_channel != filterclock_channel_name and noisesub_gate_channel != undisturbed_ok_channel and (noisesub_gate_bitmask > 0 or noisesub_gate_threshold is not None or noisesub_gate_value is not None) and noisesub_gate_channel != "None" and not noisesub_use_filter_clock:
channel_list.append((instrument, noisesub_gate_channel))
headkeys.append("noisesubgatechannel")
replace_value_list.append(0)
@@ -2652,7 +2653,7 @@ if compute_calib_statevector:
# inputs that are replaced with zeros affect h(t) for a short time before and after the zeros, so we also must account for this corrupted time.
noinvalidinput = calibration_parts.mkgate(pipeline, noinvalidinput, noinvalidinput, pow(2,no_gap_bitnum), attack_length = -int(filter_settle_time * calibstate_sr), hold_length = -int(filter_latency * calibstate_sr), queue_length = queue_length)
- #
+ #
# HW INJECTION BITS
#
@@ -2811,7 +2812,7 @@ if compute_calib_statevector:
#
# Set up gating for the power mains and noise subtraction
-if compute_calib_statevector and (any(line_witness_channel_list) or any(witness_channel_list)) and (((noisesub_gate_bitmask > 0 or noisesub_gate_threshold is not None) and noisesub_gate_channel != "None") or noisesub_use_filter_clock):
+if compute_calib_statevector and (any(line_witness_channel_list) or any(witness_channel_list)) and (((noisesub_gate_bitmask > 0 or noisesub_gate_threshold is not None or noisesub_gate_value is not None) and noisesub_gate_channel != "None") or noisesub_use_filter_clock):
if noisesub_use_filter_clock:
noisesubgate = filtersok
elif noisesub_gate_channel == obsintent_channel_name:
@@ -2830,6 +2831,10 @@ if compute_calib_statevector and (any(line_witness_channel_list) or any(witness_
noisesubgate = pipeparts.mkcapsfilter(pipeline, noisesubgate, calibstate_caps)
noisesubgate = pipeparts.mkbitvectorgen(pipeline, noisesubgate, bit_vector = pow(2,noise_sub_gate_bitnum), threshold = pow(2,filters_ok_bitnum))
noisesubgate = pipeparts.mkgeneric(pipeline, noisesubgate, "lal_logicalundersample", required_on = pow(2,noise_sub_gate_bitnum), status_out = pow(2,noise_sub_gate_bitnum), name = "noisesubgate_logicalundersample")
+ elif noisesub_gate_value is not None:
+ noisesubgate = calibration_parts.mkinsertgap(pipeline, noisesubgate, bad_data_intervals = [noisesub_gate_value, noisesub_gate_value], insert_gap = False, open_intervals = True, replace_value = 0)
+ noisesubgate = pipeparts.mkbitvectorgen(pipeline, noisesubgate, bit_vector = pow(2,noise_sub_gate_bitnum), threshold = noisesub_gate_value)
+ noisesubgate = pipeparts.mkgeneric(pipeline, noisesubgate, "lal_logicalundersample", required_on = pow(2,noise_sub_gate_bitnum), status_out = pow(2,noise_sub_gate_bitnum))
elif noisesub_gate_threshold is not None:
noisesubgate = pipeparts.mkbitvectorgen(pipeline, noisesubgate, bit_vector = pow(2,noise_sub_gate_bitnum), threshold = noisesub_gate_threshold)
noisesubgate = pipeparts.mkgeneric(pipeline, noisesubgate, "lal_logicalundersample", required_on = pow(2,noise_sub_gate_bitnum), status_out = pow(2,noise_sub_gate_bitnum))
diff --git a/gnuscripts/depcomp b/gnuscripts/depcomp
index 715e34311ed2d2dbff881aedc7e25b81db54614c..1f0aa972c93e5869704dec4354e115f498ad57bf 100755
--- a/gnuscripts/depcomp
+++ b/gnuscripts/depcomp
@@ -1,9 +1,9 @@
#! /bin/sh
# depcomp - compile a program generating dependencies as side-effects
-scriptversion=2018-03-07.03; # UTC
+scriptversion=2024-06-19.01; # UTC
-# Copyright (C) 1999-2021 Free Software Foundation, Inc.
+# Copyright (C) 1999-2024 Free Software Foundation, Inc.
# 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
@@ -47,11 +47,13 @@ Environment variables:
libtool Whether libtool is used (yes/no).
Report bugs to <bug-automake@gnu.org>.
+GNU Automake home page: <https://www.gnu.org/software/automake/>.
+General help using GNU software: <https://www.gnu.org/gethelp/>.
EOF
exit $?
;;
-v | --v*)
- echo "depcomp $scriptversion"
+ echo "depcomp (GNU Automake) $scriptversion"
exit $?
;;
esac
@@ -113,7 +115,6 @@ nl='
# These definitions help.
upper=ABCDEFGHIJKLMNOPQRSTUVWXYZ
lower=abcdefghijklmnopqrstuvwxyz
-digits=0123456789
alpha=${upper}${lower}
if test -z "$depmode" || test -z "$source" || test -z "$object"; then
@@ -128,7 +129,7 @@ tmpdepfile=${tmpdepfile-`echo "$depfile" | sed 's/\.\([^.]*\)$/.T\1/'`}
rm -f "$tmpdepfile"
-# Avoid interferences from the environment.
+# Avoid interference from the environment.
gccflag= dashmflag=
# Some modes work just like other modes, but use different flags. We
@@ -198,8 +199,8 @@ gcc3)
;;
gcc)
-## Note that this doesn't just cater to obsosete pre-3.x GCC compilers.
-## but also to in-use compilers like IMB xlc/xlC and the HP C compiler.
+## Note that this doesn't just cater to obsolete pre-3.x GCC compilers.
+## but also to in-use compilers like IBM xlc/xlC and the HP C compiler.
## (see the conditional assignment to $gccflag above).
## There are various ways to get dependency output from gcc. Here's
## why we pick this rather obscure method:
diff --git a/gst/lal/gstlal_asd.c b/gst/lal/gstlal_asd.c
index 9cc19769254cfca55a2a193754808ae9b6b155e5..92e28cd3225560fd35f568cac75d6a3aa90d09c3 100644
--- a/gst/lal/gstlal_asd.c
+++ b/gst/lal/gstlal_asd.c
@@ -253,7 +253,7 @@ static void process_input_data_ ## DTYPE(GSTLALASD *element, DTYPE *src, guint64
} \
g_free(little_asd); \
element->input_index = 0; \
- if(src_index >= element->overlap_samples) \
+ if((gint64) src_index >= element->overlap_samples) \
src_index -= element->overlap_samples; \
else { \
for(element->input_index = 0; element->input_index < element->overlap_samples - src_index; element->input_index++) \
@@ -408,7 +408,7 @@ static gboolean set_caps(GstBaseSink *sink, GstCaps *caps) {
remove(element->filename);
/* Sanity check */
- if(element->overlap_samples >= element->fft_samples) {
+ if(element->overlap_samples >= (gint64) element->fft_samples) {
GST_WARNING_OBJECT(element, "overlap_samples cannot be greater than fft_samples. Resetting overlap_samples to fft_samples / 2.");
element->overlap_samples = element->fft_samples / 2;
}
@@ -716,7 +716,7 @@ static void set_property(GObject *object, enum property id, const GValue *value,
break;
case ARG_OVERLAP_SAMPLES:
- element->overlap_samples = g_value_get_uint64(value);
+ element->overlap_samples = g_value_get_int64(value);
break;
case ARG_UPDATE_TIME:
@@ -764,7 +764,7 @@ static void get_property(GObject *object, enum property id, GValue *value, GPara
break;
case ARG_OVERLAP_SAMPLES:
- g_value_set_uint64(value, element->overlap_samples);
+ g_value_set_int64(value, element->overlap_samples);
break;
case ARG_UPDATE_TIME:
@@ -897,11 +897,11 @@ static void gstlal_asd_class_init(GSTLALASDClass *klass) {
1, G_MAXUINT64, 16384,
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
);
- properties[ARG_OVERLAP_SAMPLES] = g_param_spec_uint64(
+ properties[ARG_OVERLAP_SAMPLES] = g_param_spec_int64(
"overlap-samples",
"Overlap Samples",
"Number of input samples of overlap between consecutive FFTs",
- 0, G_MAXUINT64, 8192,
+ -G_MAXINT64, G_MAXINT64, 8192,
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
);
properties[ARG_UPDATE_TIME] = g_param_spec_uint64(
diff --git a/gst/lal/gstlal_asd.h b/gst/lal/gstlal_asd.h
index eeca4f267d3898621754a70adbfeac61a7d60bf2..e9f0324c406dbdc0095eb3d568b42cfaeae23c91 100644
--- a/gst/lal/gstlal_asd.h
+++ b/gst/lal/gstlal_asd.h
@@ -97,7 +97,7 @@ struct _GSTLALASD {
/* properties */
guint64 fft_samples;
- guint64 overlap_samples;
+ gint64 overlap_samples;
guint64 update_time;
gboolean use_td_median;
guint64 update_samples;
diff --git a/gst/lal/gstlal_insertgap.c b/gst/lal/gstlal_insertgap.c
index dbafdb5116d405d866f74787167a8537229c654a..87f914518b698aa17d032a3d93fa4d35f0fc62a2 100644
--- a/gst/lal/gstlal_insertgap.c
+++ b/gst/lal/gstlal_insertgap.c
@@ -106,7 +106,7 @@ static GstStaticPadTemplate sink_factory = GST_STATIC_PAD_TEMPLATE(
"audio/x-raw, " \
"rate = (int) [1, MAX], " \
"channels = (int) [1, MAX], " \
- "format = (string) { "GST_AUDIO_NE(U32)", " GST_AUDIO_NE(F32) ", " GST_AUDIO_NE(F64) ", " GST_AUDIO_NE(Z64) ", " GST_AUDIO_NE(Z128) " }, " \
+ "format = (string) { "GST_AUDIO_NE(S32)", "GST_AUDIO_NE(U32)", " GST_AUDIO_NE(F32) ", " GST_AUDIO_NE(F64) ", " GST_AUDIO_NE(Z64) ", " GST_AUDIO_NE(Z128) " }, " \
"layout = (string) interleaved, " \
"channel-mask = (bitmask) 0"
)
@@ -121,7 +121,7 @@ static GstStaticPadTemplate src_factory = GST_STATIC_PAD_TEMPLATE(
"audio/x-raw, " \
"rate = (int) [1, MAX], " \
"channels = (int) [1, MAX], " \
- "format = (string) { "GST_AUDIO_NE(U32)", " GST_AUDIO_NE(F32) ", " GST_AUDIO_NE(F64) ", " GST_AUDIO_NE(Z64) ", " GST_AUDIO_NE(Z128) " }, " \
+ "format = (string) { "GST_AUDIO_NE(S32)", "GST_AUDIO_NE(U32)", " GST_AUDIO_NE(F32) ", " GST_AUDIO_NE(F64) ", " GST_AUDIO_NE(Z64) ", " GST_AUDIO_NE(Z128) " }, " \
"layout = (string) interleaved, " \
"channel-mask = (bitmask) 0"
)
@@ -223,7 +223,7 @@ static GstFlowReturn push_srcbuf(GSTLALInsertGap *element, const gint8 *data, gu
#define DEFINE_CHECK_DATA(DTYPE, COMPLEX) \
-static gboolean check_data_ ## DTYPE ## COMPLEX(const DTYPE COMPLEX *data, double *bad_data_intervals, int array_length, int num_checks, gboolean remove_nan, gboolean remove_inf, gint complex_factor) { \
+static gboolean check_data_ ## DTYPE ## COMPLEX(const DTYPE COMPLEX *data, double *bad_data_intervals, int array_length, int num_checks, gboolean remove_nan, gboolean remove_inf, gint complex_factor, gboolean open_intervals) { \
int i, j; \
gboolean result_real, result_imag; \
if(complex_factor) { \
@@ -231,9 +231,16 @@ static gboolean check_data_ ## DTYPE ## COMPLEX(const DTYPE COMPLEX *data, doubl
if(bad_data_intervals) { \
result_real = FALSE; \
result_imag = FALSE; \
- for(j = 0; j < array_length; j += 4) { \
- result_real |= creal(data[i]) > bad_data_intervals[j] && creal(data[i]) < bad_data_intervals[j + 2]; \
- result_imag |= cimag(data[i]) > bad_data_intervals[j + 1] && cimag(data[i]) < bad_data_intervals[j + 3]; \
+ if(open_intervals) { \
+ for(j = 0; j < array_length; j += 4) { \
+ result_real |= creal(data[i]) >= bad_data_intervals[j] && creal(data[i]) <= bad_data_intervals[j + 2]; \
+ result_imag |= cimag(data[i]) >= bad_data_intervals[j + 1] && cimag(data[i]) <= bad_data_intervals[j + 3]; \
+ } \
+ } else { \
+ for(j = 0; j < array_length; j += 4) { \
+ result_real |= creal(data[i]) > bad_data_intervals[j] && creal(data[i]) < bad_data_intervals[j + 2]; \
+ result_imag |= cimag(data[i]) > bad_data_intervals[j + 1] && cimag(data[i]) < bad_data_intervals[j + 3]; \
+ } \
} \
if(!(result_real && result_imag)) \
return FALSE; \
@@ -247,8 +254,13 @@ static gboolean check_data_ ## DTYPE ## COMPLEX(const DTYPE COMPLEX *data, doubl
for(i = 0; i < num_checks; i++) { \
if(bad_data_intervals) { \
result_real = FALSE; \
- for(j = 0; j < array_length; j += 2) \
- result_real |= (DTYPE) data[i] > bad_data_intervals[j] && (DTYPE) data[i] < bad_data_intervals[j + 1]; \
+ if(open_intervals) { \
+ for(j = 0; j < array_length; j += 2) \
+ result_real |= (DTYPE) data[i] >= bad_data_intervals[j] && (DTYPE) data[i] <= bad_data_intervals[j + 1]; \
+ } else { \
+ for(j = 0; j < array_length; j += 2) \
+ result_real |= (DTYPE) data[i] > bad_data_intervals[j] && (DTYPE) data[i] < bad_data_intervals[j + 1]; \
+ } \
if(!result_real) \
return FALSE; \
} \
@@ -266,6 +278,7 @@ DEFINE_CHECK_DATA(float, );
DEFINE_CHECK_DATA(double, );
DEFINE_CHECK_DATA(float, complex);
DEFINE_CHECK_DATA(double, complex);
+DEFINE_CHECK_DATA(gint32, );
DEFINE_CHECK_DATA(guint32, );
@@ -343,7 +356,7 @@ static GstFlowReturn process_inbuf_ ## DTYPE ## COMPLEX(const DTYPE COMPLEX *ind
g_assert_cmpuint(max_block_length, >, 0); \
\
for(offset = 0; offset < length; offset++) { \
- data_is_bad = !check_data_ ## DTYPE ## COMPLEX(indata, element->bad_data_intervals, element->array_length, element->channels, element->remove_nan, element->remove_inf, complex_factor); \
+ data_is_bad = !check_data_ ## DTYPE ## COMPLEX(indata, element->bad_data_intervals, element->array_length, element->channels, element->remove_nan, element->remove_inf, complex_factor, element->open_intervals); \
if(data_is_bad && (element->replace_value < G_MAXDOUBLE || element->replace_value_imag < G_MAXDOUBLE)) { \
for(i = 0; i < element->channels; i++) \
outdata[i] = replace_value; \
@@ -416,11 +429,12 @@ done: \
}
-DEFINE_PROCESS_INBUF(guint32, )
-DEFINE_PROCESS_INBUF(float, )
-DEFINE_PROCESS_INBUF(double, )
-DEFINE_PROCESS_INBUF(float,complex)
-DEFINE_PROCESS_INBUF(double,complex)
+DEFINE_PROCESS_INBUF(gint32, );
+DEFINE_PROCESS_INBUF(guint32, );
+DEFINE_PROCESS_INBUF(float, );
+DEFINE_PROCESS_INBUF(double, );
+DEFINE_PROCESS_INBUF(float,complex);
+DEFINE_PROCESS_INBUF(double,complex);
static void *input_buffer_timer(void *void_element) {
@@ -480,6 +494,9 @@ naptime:
*/
GstFlowReturn result;
switch(element->data_type) {
+ case GSTLAL_INSERTGAP_S32:
+ result = process_inbuf_gint32(NULL, NULL, element, TRUE, TRUE, 0, 0, 0, ets_min, 0, FALSE);
+ break;
case GSTLAL_INSERTGAP_U32:
result = process_inbuf_guint32(NULL, NULL, element, TRUE, TRUE, 0, 0, 0, ets_min, 0, FALSE);
break;
@@ -545,7 +562,10 @@ static gboolean sink_event(GstPad *pad, GstObject *parent, GstEvent *event)
element->rate = GST_AUDIO_INFO_RATE(&info);
element->channels = GST_AUDIO_INFO_CHANNELS(&info);
element->unit_size = GST_AUDIO_INFO_BPF(&info);
- if(!strcmp(name, GST_AUDIO_NE(U32))) {
+ if(!strcmp(name, GST_AUDIO_NE(S32))) {
+ element->data_type = GSTLAL_INSERTGAP_S32;
+ g_assert_cmpuint(element->unit_size, ==, element->channels * 4);
+ } else if(!strcmp(name, GST_AUDIO_NE(U32))) {
element->data_type = GSTLAL_INSERTGAP_U32;
g_assert_cmpuint(element->unit_size, ==, element->channels * 4);
} else if(!strcmp(name, GST_AUDIO_NE(F32))) {
@@ -668,6 +688,9 @@ static GstFlowReturn chain(GstPad *pad, GstObject *parent, GstBuffer *sinkbuf)
g_mutex_lock(&element->mutex);
if(GST_BUFFER_PTS_IS_VALID(sinkbuf) && GST_BUFFER_PTS(sinkbuf) > element->last_sinkbuf_ets) {
switch(element->data_type) {
+ case GSTLAL_INSERTGAP_S32:
+ result = process_inbuf_gint32(NULL, NULL, element, TRUE, TRUE, 0, 0, 0, GST_BUFFER_PTS(sinkbuf), 0, TRUE);
+ break;
case GSTLAL_INSERTGAP_U32:
result = process_inbuf_guint32(NULL, NULL, element, TRUE, TRUE, 0, 0, 0, GST_BUFFER_PTS(sinkbuf), 0, TRUE);
break;
@@ -723,6 +746,9 @@ static GstFlowReturn chain(GstPad *pad, GstObject *parent, GstBuffer *sinkbuf)
outdata = g_malloc((sinkbuf_offset_end - sinkbuf_offset) * element->unit_size);
switch(element->data_type) {
+ case GSTLAL_INSERTGAP_S32:
+ result = process_inbuf_gint32((gint32 *) inmap.data, outdata, element, sinkbuf_gap, sinkbuf_discont, sinkbuf_offset, sinkbuf_offset_end, sinkbuf_dur, sinkbuf_pts, 0, FALSE);
+ break;
case GSTLAL_INSERTGAP_U32:
result = process_inbuf_guint32((guint32 *) inmap.data, outdata, element, sinkbuf_gap, sinkbuf_discont, sinkbuf_offset, sinkbuf_offset_end, sinkbuf_dur, sinkbuf_pts, 0, FALSE);
break;
@@ -784,6 +810,7 @@ enum property {
ARG_REPLACE_VALUE_IMAG,
ARG_KEEP_LAST_GOOD_VALUE,
ARG_BAD_DATA_INTERVALS,
+ ARG_OPEN_INTERVALS,
ARG_BLOCK_DURATION,
ARG_CHOP_LENGTH,
ARG_WAIT_TIME
@@ -838,6 +865,9 @@ static void set_property(GObject *object, enum property prop_id, const GValue *v
break;
}
+ case ARG_OPEN_INTERVALS:
+ element->open_intervals = g_value_get_boolean(value);
+ break;
case ARG_BLOCK_DURATION:
element->block_duration = g_value_get_uint64(value);
break;
@@ -901,6 +931,9 @@ static void get_property(GObject *object, enum property prop_id, GValue *value,
g_value_copy(&va, value);
g_value_unset(&va);
+ break;
+ case ARG_OPEN_INTERVALS:
+ g_value_set_boolean(value, element->open_intervals);
break;
case ARG_BLOCK_DURATION:
g_value_set_uint64(value, element->block_duration);
@@ -1091,7 +1124,7 @@ static void gstlal_insertgap_class_init(GSTLALInsertGapClass *klass)
"For a given array X[i], input data is classified \"good\" or \"bad\" as follows:\n\t\t\t"
"bad <= X[0] < good < X[1] <= bad <= X[2] < good < X[3] <= bad ...\n\t\t\t"
"For real streams, the length of this array property must be even. For complex\n\t\t\t"
- "streams, the lengthmust be a multiple of 4; even array elements (0, 2, 4, etc.)\n\t\t\t"
+ "streams, the length must be a multiple of 4; even array elements (0, 2, 4, etc.)\n\t\t\t"
"are used to check the real parts of input data and odd array elements are used\n\t\t\t"
"to check imaginary parts. If any part or channel of a sample is \"bad,\" the\n\t\t\t"
"whole sample is counted \"bad.\"",
@@ -1106,6 +1139,18 @@ static void gstlal_insertgap_class_init(GSTLALInsertGapClass *klass)
)
);
+ g_object_class_install_property(
+ gobject_class,
+ ARG_OPEN_INTERVALS,
+ g_param_spec_boolean(
+ "open-intervals",
+ "Open Intervals",
+ "If true, bad-data-intervals will be open rather than closed intervals.",
+ FALSE,
+ G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
+ )
+ );
+
g_object_class_install_property(
gobject_class,
ARG_BLOCK_DURATION,
diff --git a/gst/lal/gstlal_insertgap.h b/gst/lal/gstlal_insertgap.h
index 99e39366242ea42047163387f330fc491d109a3e..271a21ccd0d831e79d8bcd1c10d2bc6ab48d1b1e 100644
--- a/gst/lal/gstlal_insertgap.h
+++ b/gst/lal/gstlal_insertgap.h
@@ -68,7 +68,8 @@ struct _GSTLALInsertGap {
gint channels;
gint unit_size;
enum gstlal_insertgap_data_type {
- GSTLAL_INSERTGAP_U32 = 0,
+ GSTLAL_INSERTGAP_S32 = 0,
+ GSTLAL_INSERTGAP_U32,
GSTLAL_INSERTGAP_F32,
GSTLAL_INSERTGAP_F64,
GSTLAL_INSERTGAP_Z64,
@@ -97,6 +98,7 @@ struct _GSTLALInsertGap {
double replace_value_imag;
gboolean keep_last_good_value;
double *bad_data_intervals;
+ gboolean open_intervals;
gint array_length;
guint64 chop_length;
GstClockTime block_duration;
diff --git a/gst/lal/gstlal_property.c b/gst/lal/gstlal_property.c
index c404dcdcd26292723bc9298d3d822e58174785a4..a99275b6d627d037281726bc72b292ad53169751 100644
--- a/gst/lal/gstlal_property.c
+++ b/gst/lal/gstlal_property.c
@@ -89,6 +89,7 @@ enum property {
ARG_UPDATE_ABOVE_THRESHOLD,
ARG_CURRENT_AVERAGE,
ARG_TIMESTAMPED_AVERAGE,
+ ARG_UPDATE_AFTER_GAP,
ARG_FAKE
};
@@ -114,6 +115,14 @@ static void rebuild_workspace_and_reset(GObject *object) {
static void average_input_data_ ## DTYPE(GSTLALProperty *element, DTYPE *src, guint64 src_size) { \
\
gint64 i, j; \
+ if(element->update_after_gap && element->last_buffer_was_gap) { \
+ element->current_average = (double) src[0]; \
+ GST_LOG_OBJECT(element, "Just computed new property"); \
+ /* Let other elements know when the change occurred */ \
+ g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TIMESTAMPED_AVERAGE]); \
+ /* Let other elements know about the update */ \
+ g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_CURRENT_AVERAGE]); \
+ } \
if(element->update_when_change) { \
for(i = 0; i < (gint64) src_size; i++) { \
/* Check if the input value has changed */ \
@@ -135,7 +144,7 @@ static void average_input_data_ ## DTYPE(GSTLALProperty *element, DTYPE *src, gu
if(element->num_in_avg) \
start_sample = 0; \
else \
- start_sample = (gint64) (element->update_samples - (gst_util_uint64_scale_int_round(element->timestamp, element->rate, GST_SECOND) + element->average_samples - element->shift_samples) % element->update_samples) % element->update_samples; \
+ start_sample = (gint64) (element->update_samples - (gst_util_uint64_scale_int_round(element->timestamp, element->rate, GST_SECOND) + (element->average_samples - 1) - element->shift_samples) % element->update_samples) % element->update_samples; \
\
/* How many samples from this buffer will we need to add into this average? */ \
samples_to_add = element->average_samples - element->num_in_avg < (gint64) src_size - start_sample ? element->average_samples - element->num_in_avg : (gint64) src_size - start_sample; \
@@ -154,7 +163,7 @@ static void average_input_data_ ## DTYPE(GSTLALProperty *element, DTYPE *src, gu
element->current_average = element->current_sum / element->num_in_avg; \
\
/* When exactly did this change occur? */ \
- element->timestamp += gst_util_uint64_scale_int_round((guint64) (start_sample + samples_to_add), GST_SECOND, element->rate); \
+ element->timestamp += gst_util_uint64_scale_int_round((guint64) (start_sample + samples_to_add - 1), GST_SECOND, element->rate); \
\
GST_LOG_OBJECT(element, "Just computed new property"); \
\
@@ -362,11 +371,14 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
element->next_in_offset = GST_BUFFER_OFFSET_END(buffer);
GST_DEBUG_OBJECT(element, "have buffer spanning %" GST_BUFFER_BOUNDARIES_FORMAT, GST_BUFFER_BOUNDARIES_ARGS(buffer));
+ gboolean gap = GST_BUFFER_FLAG_IS_SET(buffer, GST_BUFFER_FLAG_GAP);
+
+ gst_buffer_map(buffer, &mapinfo, GST_MAP_READ);
+
/* Check if the data on this buffer is usable and if we plan to use it */
- gint64 next_start_sample = (element->update_samples - (gst_util_uint64_scale_int_round(GST_BUFFER_PTS(buffer), element->rate, GST_SECOND) + element->average_samples - element->shift_samples) % element->update_samples) % element->update_samples;
- if(!GST_BUFFER_FLAG_IS_SET(buffer, GST_BUFFER_FLAG_GAP) && mapinfo.size && (element->num_in_avg || next_start_sample < (gint64) gst_util_uint64_scale_int_round(GST_BUFFER_DURATION(buffer), element->rate, GST_SECOND) || element->update_when_change)) {
+ gint64 next_start_sample = (element->update_samples - (gst_util_uint64_scale_int_round(GST_BUFFER_PTS(buffer), element->rate, GST_SECOND) + (element->average_samples - 1) - element->shift_samples) % element->update_samples) % element->update_samples;
+ if(!gap && mapinfo.size && (element->num_in_avg || next_start_sample < (gint64) gst_util_uint64_scale_int_round(GST_BUFFER_DURATION(buffer), element->rate, GST_SECOND) || element->update_when_change)) {
/* Get the data from the buffer */
- gst_buffer_map(buffer, &mapinfo, GST_MAP_READ);
switch(element->data_type) {
case GSTLAL_PROPERTY_U8:
@@ -397,8 +409,10 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
g_assert_not_reached();
break;
}
- gst_buffer_unmap(buffer, &mapinfo);
}
+ gst_buffer_unmap(buffer, &mapinfo);
+
+ element->last_buffer_was_gap = gap;
return result;
}
@@ -453,6 +467,10 @@ static void set_property(GObject *object, enum property id, const GValue *value,
element->current_average = g_value_get_double(value);
break;
+ case ARG_UPDATE_AFTER_GAP:
+ element->update_after_gap = g_value_get_boolean(value);
+ break;
+
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID(object, id, pspec);
break;
@@ -514,6 +532,10 @@ static void get_property(GObject *object, enum property id, GValue *value, GPara
g_value_unset(&varray);
break;
+ case ARG_UPDATE_AFTER_GAP:
+ g_value_set_boolean(value, element->update_after_gap);
+ break;
+
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID(object, id, pspec);
break;
@@ -630,7 +652,7 @@ static void gstlal_property_class_init(GSTLALPropertyClass *klass) {
"average-samples",
"Average Samples",
"Number of input samples to average before updating the property",
- 0, G_MAXINT64, 1,
+ 1, G_MAXINT64, 1,
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
);
properties[ARG_SHIFT_SAMPLES] = g_param_spec_int64(
@@ -684,6 +706,13 @@ static void gstlal_property_class_init(GSTLALPropertyClass *klass) {
),
G_PARAM_READABLE | G_PARAM_STATIC_STRINGS
);
+ properties[ARG_UPDATE_AFTER_GAP] = g_param_spec_boolean(
+ "update-after-gap",
+ "Update After Gap",
+ "If true, updates will happen anytime a gap in the data ends.",
+ FALSE,
+ G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
+ );
g_object_class_install_property(
@@ -726,6 +755,11 @@ static void gstlal_property_class_init(GSTLALPropertyClass *klass) {
ARG_TIMESTAMPED_AVERAGE,
properties[ARG_TIMESTAMPED_AVERAGE]
);
+ g_object_class_install_property(
+ gobject_class,
+ ARG_UPDATE_AFTER_GAP,
+ properties[ARG_UPDATE_AFTER_GAP]
+ );
}
diff --git a/gst/lal/gstlal_property.h b/gst/lal/gstlal_property.h
index c18c2a0dd4756bb628d25191026ebc2285895bf8..110bc3d3d8e4dc3d06a57bfa221c6034254cb4a4 100644
--- a/gst/lal/gstlal_property.h
+++ b/gst/lal/gstlal_property.h
@@ -75,6 +75,7 @@ struct _GSTLALProperty {
/* filter memory */
gint64 num_in_avg;
double current_sum;
+ gboolean last_buffer_was_gap;
union {
struct {
@@ -112,6 +113,7 @@ struct _GSTLALProperty {
gboolean update_above_threshold;
double current_average;
guint64 timestamp;
+ gboolean update_after_gap;
};
diff --git a/gst/lal/gstlal_smoothkappas.c b/gst/lal/gstlal_smoothkappas.c
index a62a59b8a4eea5e24a4742f6d40b8c20ed091d50..6cedb348618d0860eb27c67c51c651237104f16b 100644
--- a/gst/lal/gstlal_smoothkappas.c
+++ b/gst/lal/gstlal_smoothkappas.c
@@ -151,8 +151,8 @@ static void set_metadata(GSTLALSmoothKappas *element, GstBuffer *buf, guint64 ou
GST_BUFFER_OFFSET(buf) = element->next_out_offset;
element->next_out_offset += outsamples;
GST_BUFFER_OFFSET_END(buf) = element->next_out_offset;
- GST_BUFFER_TIMESTAMP(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET(buf) - element->offset0, GST_SECOND, element->rate);
- GST_BUFFER_DURATION(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET_END(buf) - element->offset0, GST_SECOND, element->rate) - GST_BUFFER_TIMESTAMP(buf);
+ GST_BUFFER_PTS(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET(buf) - element->offset0, GST_SECOND, element->rate);
+ GST_BUFFER_DURATION(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET_END(buf) - element->offset0, GST_SECOND, element->rate) - GST_BUFFER_PTS(buf);
GST_BUFFER_FLAG_UNSET(buf, GST_BUFFER_FLAG_GAP);
if(G_UNLIKELY(element->need_discont)) {
GST_BUFFER_FLAG_SET(buf, GST_BUFFER_FLAG_DISCONT);
@@ -213,12 +213,16 @@ static void get_new_median(double new_element, double *fifo_array, double *curre
else if(number_less > num_samples / 2)
*current_median = less;
} else {
- if(number_greater >= 1 + num_samples / 2)
+ if(number_greater > num_samples / 2)
*current_median = (Greater + greater) / 2;
- else if(number_less >= 1 + num_samples / 2)
+ else if(number_less > num_samples / 2)
*current_median = (Less + less) / 2;
else if(number_greater == num_samples / 2 && number_less == num_samples / 2)
*current_median = (greater + less) / 2;
+ else if(number_greater == num_samples / 2)
+ *current_median = (greater + *current_median) / 2;
+ else if(number_less == num_samples / 2)
+ *current_median = (less + *current_median) / 2;
}
return;
@@ -238,7 +242,7 @@ static double get_average(double new_element, double *fifo_array, gint array_siz
double sum = 0;
int i, ii, i_start, num_samples;
num_samples = (no_default && samples_in_filter < array_size) ? samples_in_filter : array_size;
- i_start = array_is_imaginary ? *index_im : *index_re + array_size - num_samples;
+ i_start = (array_is_imaginary ? *index_im : *index_re) + array_size - num_samples;
for(i = i_start; i < i_start + num_samples; i++) {
ii = i % array_size;
sum += fifo_array[ii];
@@ -247,6 +251,121 @@ static double get_average(double new_element, double *fifo_array, gint array_siz
}
+static void reset(GSTLALSmoothKappas *element) {
+ element->samples_in_filter = element->min_array_size;
+ int i, num_default, ii, i_start, i_stop, num_less, num_greater, num_equal;
+ double greater, less, replace_value_re, replace_value_im;
+
+ /* Update median array and bookkeeping */
+ if(element->array_size > element->samples_in_filter) {
+ num_default = element->array_size - element->min_array_size;
+ /* Replace removed values with default values */
+ if(!element->default_to_median) {
+ for(i = element->index_re; i < element->index_re + num_default; i++)
+ element->fifo_array_re[i % element->array_size] = element->default_kappa_re;
+ for(i = element->index_im; i < element->index_im + num_default; i++)
+ element->fifo_array_im[i % element->array_size] = element->default_kappa_im;
+ }
+ /* Find new median */
+ num_less = 0;
+ num_equal = 0;
+ num_greater = 0;
+ greater = G_MAXDOUBLE;
+ less = -G_MAXDOUBLE;
+ /* Real part */
+ while(num_less + num_equal + num_greater == 0 || num_greater > element->array_size / 2 || num_less > element->array_size / 2) {
+ num_less = 0;
+ num_equal = 0;
+ num_greater = 0;
+ i_start = (element->index_re + element->array_size - element->samples_in_filter) % element->array_size;
+ i_stop = i_start + (element->no_default ? element->samples_in_filter : element->array_size);
+ for(i = i_start; i < i_stop; i++) {
+ ii = i % element->array_size;
+ if(element->fifo_array_re[ii] < element->current_median_re) {
+ num_less++;
+ if(element->fifo_array_re[ii] > less)
+ less = element->fifo_array_re[ii];
+ } else if(element->fifo_array_re[ii] == element->current_median_re)
+ num_equal++;
+ else {
+ num_greater++;
+ if(element->fifo_array_re[ii] < greater)
+ greater = element->fifo_array_re[ii];
+ }
+ }
+ if(num_less == element->array_size / 2 && num_greater == element->array_size / 2 && !(element->array_size % 2))
+ element->current_median_re = (less + greater) / 2;
+ else if(num_less == element->array_size / 2 && !(element->array_size % 2))
+ element->current_median_re = (less + element->current_median_re) / 2;
+ else if(num_greater == element->array_size / 2 && !(element->array_size % 2))
+ element->current_median_re = (element->current_median_re + greater) / 2;
+ else if(num_greater > element->array_size / 2)
+ element->current_median_re = greater;
+ else if(num_less > element->array_size / 2)
+ element->current_median_re = less;
+ /* In all other cases, do nothing */
+ }
+ num_less = 0;
+ num_equal = 0;
+ num_greater = 0;
+ greater = G_MAXDOUBLE;
+ less = -G_MAXDOUBLE;
+ /* Imaginary part */
+ while(num_less + num_equal + num_greater == 0 || num_greater > element->array_size / 2 || num_less > element->array_size / 2) {
+ num_less = 0;
+ num_equal = 0;
+ num_greater = 0;
+ i_start = (element->index_im + element->array_size - element->samples_in_filter) % element->array_size;
+ i_stop = i_start + (element->no_default ? element->samples_in_filter : element->array_size);
+ for(i = i_start; i < i_stop; i++) {
+ ii = i % element->array_size;
+ if(element->fifo_array_im[ii] < element->current_median_im) {
+ num_less++;
+ if(element->fifo_array_im[ii] > less)
+ less = element->fifo_array_im[ii];
+ } else if(element->fifo_array_im[ii] == element->current_median_im)
+ num_equal++;
+ else {
+ num_greater++;
+ if(element->fifo_array_im[ii] < greater)
+ greater = element->fifo_array_im[ii];
+ }
+ }
+ if(num_less == element->array_size / 2 && num_greater == element->array_size / 2 && !(element->array_size % 2))
+ element->current_median_im = (less + greater) / 2;
+ else if(num_less == element->array_size / 2 && !(element->array_size % 2))
+ element->current_median_im = (less + element->current_median_im) / 2;
+ else if(num_greater == element->array_size / 2 && !(element->array_size % 2))
+ element->current_median_im = (element->current_median_im + greater) / 2;
+ else if(num_greater > element->array_size / 2)
+ element->current_median_im = greater;
+ else if(num_less > element->array_size / 2)
+ element->current_median_im = less;
+ /* In all other cases, do nothing */
+ }
+ if(element->default_to_median) {
+ /* Fill removed values with current median */
+ for(i = element->index_re; i < element->index_re + num_default; i++)
+ element->fifo_array_re[i % element->array_size] = element->current_median_re;
+ for(i = element->index_im; i < element->index_im + num_default; i++)
+ element->fifo_array_im[i % element->array_size] = element->current_median_im;
+ }
+ }
+ /* Update average array and bookkeeping */
+ if(element->avg_array_size > element->samples_in_filter) {
+ num_default = element->avg_array_size - element->min_array_size;
+ replace_value_re = element->default_to_median ? element->current_median_re : element->default_kappa_re;
+ replace_value_im = element->default_to_median ? element->current_median_im : element->default_kappa_im;
+ for(i = element->avg_index_re; i < element->avg_index_re + num_default; i++)
+ element->avg_array_re[i % element->avg_array_size] = replace_value_re;
+ for(i = element->avg_index_im; i < element->avg_index_im + num_default; i++)
+ element->avg_array_im[i % element->avg_array_size] = replace_value_im;
+ }
+
+ return;
+}
+
+
#define DEFINE_SMOOTH_BUFFER(DTYPE) \
static GstFlowReturn smooth_buffer_ ## DTYPE(const DTYPE *src, guint64 src_size, DTYPE *dst, guint64 dst_size, double *fifo_array, double *avg_array, double default_kappa, double *current_median, double maximum_offset, gint array_size, gint avg_array_size, int *index_re, int *index_im, int *avg_index_re, int *avg_index_im, int *num_bad_in_avg_re, gboolean gap, gboolean default_to_median, gboolean track_bad_kappa, int *samples_in_filter, gboolean no_default, guint64 *total_samples, gboolean reject_zeros) { \
guint64 i; \
@@ -614,9 +733,10 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
GSTLALSmoothKappas *element = GSTLAL_SMOOTHKAPPAS(trans);
GstMapInfo inmap, outmap;
GstFlowReturn result;
+ int i, waste_samples;
/* How many samples do we need to throw away based on the filter latency? */
- int waste_samples = (int) (element->filter_latency * (element->array_size + element->avg_array_size - 2));
+ waste_samples = (int) (element->filter_latency * (element->array_size + element->avg_array_size - 2));
/*
* Check for discontinuity
@@ -639,7 +759,6 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
element->current_median_re = element->default_kappa_re;
element->current_median_im = element->default_kappa_im;
- int i;
for(i = 0; i < element->array_size; i++) {
element->fifo_array_re[i] = element->default_kappa_re;
element->fifo_array_im[i] = element->default_kappa_im;
@@ -649,6 +768,19 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
element->avg_array_im[i] = element->default_kappa_im;
}
}
+ /* Reset the running median and mean calculation if requested at this time */
+ g_mutex_lock(&element->reset_mutex);
+ if(element->num_reset_times > 0) {
+ if(GST_BUFFER_PTS(inbuf) >= *element->reset_times) {
+ if(element->min_array_size < element->samples_in_filter)
+ reset(element);
+ element->num_reset_times--;
+ for(i = 0; (guint) i < element->num_reset_times; i++)
+ element->reset_times[i] = element->reset_times[i + 1];
+ element->reset_times = g_realloc(element->reset_times, element->num_reset_times * sizeof(guint64));
+ }
+ }
+ g_mutex_unlock(&element->reset_mutex);
gst_buffer_map(outbuf, &outmap, GST_MAP_WRITE);
gst_buffer_map(inbuf, &inmap, GST_MAP_READ);
@@ -711,7 +843,7 @@ enum property {
ARG_FILTER_LATENCY,
ARG_NO_DEFAULT,
ARG_CLEAR_ON_GAP,
- ARG_RESET,
+ ARG_RESET_TIME,
ARG_REJECT_ZEROS
};
@@ -759,120 +891,32 @@ static void set_property(GObject *object, enum property prop_id, const GValue *v
case ARG_CLEAR_ON_GAP:
element->clear_on_gap = g_value_get_boolean(value);
break;
- case ARG_RESET:
- element->reset = g_value_get_boolean(value);
- if(element->reset) {
- if(element->min_array_size < element->samples_in_filter) {
- element->samples_in_filter = element->min_array_size;
- /* Update median array and bookkeeping */
- if(element->array_size > element->samples_in_filter) {
- int i, num_default = element->array_size - element->min_array_size;
- /* Replace removed values with default values */
- if(!element->default_to_median) {
- for(i = element->index_re; i < element->index_re + num_default; i++)
- element->fifo_array_re[i % element->array_size] = element->default_kappa_re;
- for(i = element->index_im; i < element->index_im + num_default; i++)
- element->fifo_array_im[i % element->array_size] = element->default_kappa_im;
- }
- /* Find new median */
- int ii, i_start, i_stop, num_less, num_greater, num_equal;
- double greater, less;
- num_less = 0;
- num_equal = 0;
- num_greater = 0;
- greater = G_MAXDOUBLE;
- less = -G_MAXDOUBLE;
- /* Real part */
- while(num_less + num_equal + num_greater == 0 || num_greater > element->array_size / 2 || num_less > element->array_size / 2) {
- num_less = 0;
- num_equal = 0;
- num_greater = 0;
- i_start = (element->index_re + element->array_size - element->samples_in_filter) % element->array_size;
- i_stop = i_start + (element->no_default ? element->samples_in_filter : element->array_size);
- for(i = i_start; i < i_stop; i++) {
- ii = i % element->array_size;
- if(element->fifo_array_re[ii] < element->current_median_re) {
- num_less++;
- if(element->fifo_array_re[ii] > less)
- less = element->fifo_array_re[ii];
- } else if(element->fifo_array_re[ii] == element->current_median_re)
- num_equal++;
- else {
- num_greater++;
- if(element->fifo_array_re[ii] < greater)
- greater = element->fifo_array_re[ii];
- }
- }
- if(num_less == element->array_size / 2 && num_greater == element->array_size / 2 && !(element->array_size % 2))
- element->current_median_re = (less + greater) / 2;
- else if(num_less == element->array_size / 2 && !(element->array_size % 2))
- element->current_median_re = (less + element->current_median_re) / 2;
- else if(num_greater == element->array_size / 2 && !(element->array_size % 2))
- element->current_median_re = (element->current_median_re + greater) / 2;
- else if(num_greater > element->array_size / 2)
- element->current_median_re = greater;
- else if(num_less > element->array_size / 2)
- element->current_median_re = less;
- /* In all other cases, do nothing */
- }
- num_less = 0;
- num_equal = 0;
- num_greater = 0;
- greater = G_MAXDOUBLE;
- less = -G_MAXDOUBLE;
- /* Imaginary part */
- while(num_less + num_equal + num_greater == 0 || num_greater > element->array_size / 2 || num_less > element->array_size / 2) {
- num_less = 0;
- num_equal = 0;
- num_greater = 0;
- i_start = (element->index_im + element->array_size - element->samples_in_filter) % element->array_size;
- i_stop = i_start + (element->no_default ? element->samples_in_filter : element->array_size);
- for(i = i_start; i < i_stop; i++) {
- ii = i % element->array_size;
- if(element->fifo_array_im[ii] < element->current_median_im) {
- num_less++;
- if(element->fifo_array_im[ii] > less)
- less = element->fifo_array_im[ii];
- } else if(element->fifo_array_im[ii] == element->current_median_im)
- num_equal++;
- else {
- num_greater++;
- if(element->fifo_array_im[ii] < greater)
- greater = element->fifo_array_im[ii];
- }
- }
- if(num_less == element->array_size / 2 && num_greater == element->array_size / 2 && !(element->array_size % 2))
- element->current_median_im = (less + greater) / 2;
- else if(num_less == element->array_size / 2 && !(element->array_size % 2))
- element->current_median_im = (less + element->current_median_im) / 2;
- else if(num_greater == element->array_size / 2 && !(element->array_size % 2))
- element->current_median_im = (element->current_median_im + greater) / 2;
- else if(num_greater > element->array_size / 2)
- element->current_median_im = greater;
- else if(num_less > element->array_size / 2)
- element->current_median_im = less;
- /* In all other cases, do nothing */
- }
- if(element->default_to_median) {
- /* Fill removed values with current median */
- for(i = element->index_re; i < element->index_re + num_default; i++)
- element->fifo_array_re[i % element->array_size] = element->current_median_re;
- for(i = element->index_im; i < element->index_im + num_default; i++)
- element->fifo_array_im[i % element->array_size] = element->current_median_im;
- }
- }
- /* Update average array and bookkeeping */
- if(element->avg_array_size > element->samples_in_filter) {
- int i, num_default = element->avg_array_size - element->min_array_size;
- double replace_value_re = element->default_to_median ? element->current_median_re : element->default_kappa_re;
- double replace_value_im = element->default_to_median ? element->current_median_im : element->default_kappa_im;
- for(i = element->avg_index_re; i < element->avg_index_re + num_default; i++)
- element->avg_array_re[i % element->avg_array_size] = replace_value_re;
- for(i = element->avg_index_im; i < element->avg_index_im + num_default; i++)
- element->avg_array_im[i % element->avg_array_size] = replace_value_im;
- }
+ case ARG_RESET_TIME: ;
+ guint64 old_reset_time, new_reset_time = g_value_get_uint64(value);
+ guint i = 0;
+ /* Only store nonzero reset times that we don't already have */
+ gboolean store_new_reset_time = new_reset_time > 0;
+ if(store_new_reset_time) {
+ for(i = 0; i < element->num_reset_times; i++)
+ store_new_reset_time = store_new_reset_time && new_reset_time != element->reset_times[i];
+ }
+ if(store_new_reset_time) {
+ g_mutex_lock(&element->reset_mutex);
+ element->num_reset_times++;
+ element->reset_times = g_realloc(element->reset_times, element->num_reset_times * sizeof(guint64));
+ /* Store the reset times in increasing order */
+ while(element->reset_times[i] < new_reset_time && i < element->num_reset_times - 1)
+ i++;
+ old_reset_time = element->reset_times[i];
+ element->reset_times[i] = new_reset_time;
+ /* Shift values in the array if we had to insert the new reset time in the middle */
+ while(i < element->num_reset_times - 1) {
+ i++;
+ new_reset_time = element->reset_times[i];
+ element->reset_times[i] = old_reset_time;
+ old_reset_time = new_reset_time;
}
- element->reset = FALSE;
+ g_mutex_unlock(&element->reset_mutex);
}
break;
case ARG_REJECT_ZEROS:
@@ -930,8 +974,11 @@ static void get_property(GObject *object, enum property prop_id, GValue *value,
case ARG_CLEAR_ON_GAP:
g_value_set_boolean(value, element->clear_on_gap);
break;
- case ARG_RESET:
- g_value_set_boolean(value, element->reset);
+ case ARG_RESET_TIME:
+ if(element->num_reset_times > 0)
+ g_value_set_uint64(value, element->reset_times[0]);
+ else
+ g_value_set_uint64(value, 0);
break;
case ARG_REJECT_ZEROS:
g_value_set_boolean(value, element->reject_zeros);
@@ -963,6 +1010,12 @@ static void finalize(GObject *object)
element->avg_array_re = NULL;
g_free(element->avg_array_im);
element->avg_array_im = NULL;
+ if(element->reset_times) {
+ g_free(element->reset_times);
+ element->reset_times = NULL;
+ element->num_reset_times = 0;
+ }
+ g_mutex_clear(&element->reset_mutex);
G_OBJECT_CLASS(gstlal_smoothkappas_parent_class)->finalize(object);
}
@@ -1154,13 +1207,14 @@ static void gstlal_smoothkappas_class_init(GSTLALSmoothKappasClass *klass)
);
g_object_class_install_property(
gobject_class,
- ARG_RESET,
- g_param_spec_boolean(
- "reset",
- "Reset",
- "If set to true, values present in the median and average arrays will be\n\t\t\t"
- "deleted (min array sizes are kept), and the calculation will start over.",
- FALSE,
+ ARG_RESET_TIME,
+ g_param_spec_uint64(
+ "reset-time",
+ "Reset Time",
+ "Next presentation timestamp (in ns) at which values present in the median\n\t\t\t"
+ "and average arrays will be deleted (min array sizes are kept), and the\n\t\t\t"
+ "calculation will start over. Default is to not reset.",
+ 0, G_MAXUINT64, 0,
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT | GST_PARAM_CONTROLLABLE
)
);
@@ -1209,6 +1263,11 @@ static void gstlal_smoothkappas_init(GSTLALSmoothKappas *element) {
element->num_bad_in_avg_im = G_MAXINT;
element->samples_in_filter = 0;
element->total_samples = 0;
+ element->reset_times = NULL;
+ element->num_reset_times = 0;
+
+ g_mutex_init(&element->reset_mutex);
+
gst_base_transform_set_qos_enabled(GST_BASE_TRANSFORM(element), TRUE);
gst_base_transform_set_gap_aware(GST_BASE_TRANSFORM(element), TRUE);
}
diff --git a/gst/lal/gstlal_smoothkappas.h b/gst/lal/gstlal_smoothkappas.h
index 0fba4d468209c5d4481836e6d55495c322c815bb..434664d24342078a56f74bb34902a2da531f3b18 100644
--- a/gst/lal/gstlal_smoothkappas.h
+++ b/gst/lal/gstlal_smoothkappas.h
@@ -96,6 +96,8 @@ struct _GSTLALSmoothKappas {
int samples_in_filter;
guint64 total_samples;
+ GMutex reset_mutex;
+
/* properties */
int array_size;
int min_array_size;
@@ -109,7 +111,8 @@ struct _GSTLALSmoothKappas {
double filter_latency;
gboolean no_default;
gboolean clear_on_gap;
- gboolean reset;
+ guint64 *reset_times;
+ guint num_reset_times;
gboolean reject_zeros;
};
diff --git a/gst/lal/gstlal_stdev.c b/gst/lal/gstlal_stdev.c
index b6482d1b4a1258039292349a2e4b4150f65d85c3..69ef5c34de28998a500030408323e09d0ab30608 100644
--- a/gst/lal/gstlal_stdev.c
+++ b/gst/lal/gstlal_stdev.c
@@ -883,6 +883,23 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
g_assert_cmpuint(inmap.size % element->unit_size_in, ==, 0);
g_assert_cmpuint(outmap.size % element->unit_size_out, ==, 0);
+ /* Reset the standard deviation calculation if requested at this time */
+ g_mutex_lock(&element->reset_mutex);
+ if(element->num_reset_times > 0) {
+ if(GST_BUFFER_PTS(inbuf) >= *element->reset_times) {
+ if(element->min_array_size < element->samples_in_array) {
+ element->samples_in_array = element->min_array_size;
+ element->start_index = (element->array_size + element->array_index - element->samples_in_array) % element->array_size;
+ }
+ element->num_reset_times--;
+ guint i;
+ for(i = 0; i < element->num_reset_times; i++)
+ element->reset_times[i] = element->reset_times[i + 1];
+ element->reset_times = g_realloc(element->reset_times, element->num_reset_times * sizeof(guint64));
+ }
+ }
+ g_mutex_unlock(&element->reset_mutex);
+
/* Process data in buffer */
switch(element->data_type) {
case GSTLAL_STDEV_F32:
@@ -937,7 +954,7 @@ enum property {
ARG_PERCENTILE,
ARG_FILTER_LATENCY,
ARG_MIN_ARRAY_SIZE,
- ARG_RESET
+ ARG_RESET_TIME
};
@@ -966,14 +983,32 @@ static void set_property(GObject *object, enum property prop_id, const GValue *v
case ARG_MIN_ARRAY_SIZE:
element->min_array_size = g_value_get_uint64(value);
break;
- case ARG_RESET:
- element->reset = g_value_get_boolean(value);
- if(element->reset) {
- if(element->min_array_size < element->samples_in_array) {
- element->samples_in_array = element->min_array_size;
- element->start_index = (element->array_size + element->array_index - element->samples_in_array) % element->array_size;
+ case ARG_RESET_TIME: ;
+ guint64 old_reset_time, new_reset_time = g_value_get_uint64(value);
+ guint i = 0;
+ /* Only store nonzero reset times that we don't already have */
+ gboolean store_new_reset_time = new_reset_time > 0;
+ if(store_new_reset_time) {
+ for(i = 0; i < element->num_reset_times; i++)
+ store_new_reset_time = store_new_reset_time && new_reset_time != element->reset_times[i];
+ }
+ if(store_new_reset_time) {
+ g_mutex_lock(&element->reset_mutex);
+ element->num_reset_times++;
+ element->reset_times = g_realloc(element->reset_times, element->num_reset_times * sizeof(guint64));
+ /* Store the reset times in increasing order */
+ while(element->reset_times[i] < new_reset_time && i < element->num_reset_times - 1)
+ i++;
+ old_reset_time = element->reset_times[i];
+ element->reset_times[i] = new_reset_time;
+ /* Shift values in the array if we had to insert the new reset time in the middle */
+ while(i < element->num_reset_times - 1) {
+ i++;
+ new_reset_time = element->reset_times[i];
+ element->reset_times[i] = old_reset_time;
+ old_reset_time = new_reset_time;
}
- element->reset = FALSE;
+ g_mutex_unlock(&element->reset_mutex);
}
break;
default:
@@ -1010,8 +1045,11 @@ static void get_property(GObject *object, enum property prop_id, GValue *value,
case ARG_MIN_ARRAY_SIZE:
g_value_set_uint64(value, element->min_array_size);
break;
- case ARG_RESET:
- g_value_set_boolean(value, element->reset);
+ case ARG_RESET_TIME:
+ if(element->num_reset_times > 0)
+ g_value_set_uint64(value, element->reset_times[0]);
+ else
+ g_value_set_uint64(value, 0);
break;
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec);
@@ -1033,6 +1071,12 @@ static void finalize(GObject *object) {
free_workspace(element);
gst_object_unref(element->srcpad);
element->srcpad = NULL;
+ if(element->reset_times) {
+ g_free(element->reset_times);
+ element->reset_times = NULL;
+ element->num_reset_times = 0;
+ }
+ g_mutex_clear(&element->reset_mutex);
G_OBJECT_CLASS(gstlal_stdev_parent_class)->finalize(object);
}
@@ -1148,13 +1192,14 @@ static void gstlal_stdev_class_init(GSTLALStDevClass *klass) {
);
g_object_class_install_property(
gobject_class,
- ARG_RESET,
- g_param_spec_boolean(
- "reset",
- "Reset",
- "If set to true, values present in the standard deviation array will be\n\t\t\t"
- "deleted, and the calculation will start over.",
- FALSE,
+ ARG_RESET_TIME,
+ g_param_spec_uint64(
+ "reset-time",
+ "Reset Time",
+ "Next presentation timestamp (in ns) at which values present in the standard\n\t\t\t"
+ "deviation array will be deleted (min array sizes are kept), and the\n\t\t\t"
+ "calculation will start over. Default is to not reset.",
+ 0, G_MAXUINT64, 0,
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT | GST_PARAM_CONTROLLABLE
)
);
@@ -1184,6 +1229,11 @@ static void gstlal_stdev_init(GSTLALStDev *element) {
element->samples_in_array = 0;
element->total_insamples = 0;
element->need_discont = TRUE;
+ element->reset_times = NULL;
+ element->num_reset_times = 0;
+
+ g_mutex_init(&element->reset_mutex);
+
gst_base_transform_set_qos_enabled(GST_BASE_TRANSFORM(element), TRUE);
gst_base_transform_set_gap_aware(GST_BASE_TRANSFORM(element), TRUE);
}
diff --git a/gst/lal/gstlal_stdev.h b/gst/lal/gstlal_stdev.h
index 07aed840b451b1d92e4fff4d2f25411a304df2ac..090e1b44150171e2f74e9f1e53c6fac5ebaab847 100644
--- a/gst/lal/gstlal_stdev.h
+++ b/gst/lal/gstlal_stdev.h
@@ -131,6 +131,8 @@ struct _GSTLALStDev {
guint64 buffer_index;
guint64 samples_in_array;
+ GMutex reset_mutex;
+
/* properties */
guint64 array_size;
guint64 coherence_length;
@@ -138,7 +140,8 @@ struct _GSTLALStDev {
double percentile;
double filter_latency;
guint64 min_array_size;
- gboolean reset;
+ guint64 *reset_times;
+ guint num_reset_times;
};
diff --git a/lib/gstlal_firtools.c b/lib/gstlal_firtools.c
index 42265e20648621c9a5e9363d15797ec140a7c7cf..062f44b832b884702e3503cdb7d8873cf6675472 100644
--- a/lib/gstlal_firtools.c
+++ b/lib/gstlal_firtools.c
@@ -2727,14 +2727,14 @@ FREQRESP(long, double);
#define GSTLAL_ASD(LONG, DTYPE, LLL, FF) \
-LONG DTYPE *gstlal_asd_ ## LONG ## DTYPE(LONG DTYPE *data, guint64 N, gint rate, guint64 fft_samples, guint64 fft_overlap, LONG DTYPE *window) { \
+LONG DTYPE *gstlal_asd_ ## LONG ## DTYPE(LONG DTYPE *data, guint64 N, gint rate, guint64 fft_samples, gint64 fft_overlap, LONG DTYPE *window) { \
\
guint64 num_ffts, i, j, j_start, asd_samples; \
LONG DTYPE *fft_input, *asd, win_sum; \
complex LONG DTYPE *fft; \
\
/* Sanity check */ \
- if(fft_overlap >= fft_samples) \
+ if(fft_overlap >= (gint64) fft_samples) \
return NULL; \
\
/* How many FFTs will we take? */ \
diff --git a/lib/gstlal_firtools.h b/lib/gstlal_firtools.h
index aca03839fadd55a67c05cd270323ff834f15d73d..f35b8cc3af3294dba78edc3aab6c562c81b7117f 100644
--- a/lib/gstlal_firtools.h
+++ b/lib/gstlal_firtools.h
@@ -330,11 +330,11 @@ complex double *freqresp_double(double *filt, guint N, guint delay_samples, guin
long complex double *freqresp_longdouble(long double *filt, guint N, guint delay_samples, guint samples_per_lobe);
-float *gstlal_asd_float(float *data, guint64 N, gint rate, guint64 fft_samples, guint64 fft_overlap, float *window);
+float *gstlal_asd_float(float *data, guint64 N, gint rate, guint64 fft_samples, gint64 fft_overlap, float *window);
-double *gstlal_asd_double(double *data, guint64 N, gint rate, guint64 fft_samples, guint64 fft_overlap, double *window);
+double *gstlal_asd_double(double *data, guint64 N, gint rate, guint64 fft_samples, gint64 fft_overlap, double *window);
-long double *gstlal_asd_longdouble(long double *data, guint64 N, gint rate, guint64 fft_samples, guint64 fft_overlap, long double *window);
+long double *gstlal_asd_longdouble(long double *data, guint64 N, gint rate, guint64 fft_samples, gint64 fft_overlap, long double *window);
G_END_DECLS
diff --git a/python/calibration_parts.py b/python/calibration_parts.py
index 602501e6b54c92768eceb6b2e87e595477f7e84c..c635c29d14de7e9149c90ce9ff18ed1ae34f58b8 100644
--- a/python/calibration_parts.py
+++ b/python/calibration_parts.py
@@ -375,7 +375,7 @@ def hook_up(pipeline, demux, channel_name, instrument, buffer_length, element_na
gap = True
replace_value = 0
- head = mkinsertgap(pipeline, None, bad_data_intervals = [-input_max, -input_min, input_min, input_max], insert_gap = gap, remove_gap = not gap, fill_discont = True, block_duration = int(1000000000 * buffer_length), replace_value = replace_value, name = "insertgap_%s%s" % (channel_name, element_name_suffix), wait_time = int(1000000000 * wait_time), keep_last_good_value = keep_last_good_value)
+ head = mkinsertgap(pipeline, None, bad_data_intervals = [-input_max, -input_min, 0, 0, input_min, input_max], open_intervals = True, insert_gap = gap, remove_gap = not gap, fill_discont = True, block_duration = int(1000000000 * buffer_length), replace_value = replace_value, name = "insertgap_%s%s" % (channel_name, element_name_suffix), wait_time = int(1000000000 * wait_time), keep_last_good_value = keep_last_good_value)
pipeparts.src_deferred_link(demux, "%s:%s" % (instrument, channel_name), head.get_static_pad("sink"))
# The following queue is needed to store data that sometimes accumulates in each channel
# at stream start when some elements are not yet ready to process data. Many of these
@@ -583,12 +583,15 @@ def remove_lines_with_witnesses(pipeline, signal, witnesses, freqs, freq_vars, f
ones = pipeparts.mktee(pipeline, mkpow(pipeline, line_in_signal, exponent = 0.0))
line_in_witnesses = []
+ witness_downsampled = []
witness_exists = []
tfs_at_f = [None] * len(witnesses[m]) * (len(witnesses[m]) + 1)
for i in range(len(witnesses[m])):
# Check whether witness exists and use to gate TF
- wit_exists = mkresample(pipeline, witnesses[m][i], 0, False, compute_rate)
- wit_exists = mkinsertgap(pipeline, wit_exists, bad_data_intervals = [-input_max, 0.9999 * witness_absent_value, 1.0001 * witness_absent_value, input_max])
+ wit_downsampled = mkresample(pipeline, witnesses[m][i], 0, False, compute_rate)
+ wit_downsampled = pipeparts.mktee(pipeline, wit_downsampled)
+ witness_downsampled.append(wit_downsampled)
+ wit_exists = mkinsertgap(pipeline, wit_downsampled, bad_data_intervals = [-input_max, 0.9999 * witness_absent_value, 1.0001 * witness_absent_value, input_max])
wit_exists = pipeparts.mkbitvectorgen(pipeline, wit_exists, nongap_is_control = True, bit_vector = 1)
wit_exists = pipeparts.mktee(pipeline, pipeparts.mkcapsfilter(pipeline, wit_exists, "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % compute_rate))
witness_exists.append(wit_exists)
@@ -613,10 +616,11 @@ def remove_lines_with_witnesses(pipeline, signal, witnesses, freqs, freq_vars, f
# gst_audio_format_info_fill_silence: assertion 'dest != NULL' failed"
# if a GAP flag is set on a buffer received by either sink pad.
tf_at_f = complex_division(pipeline, line_in_signal, line_in_witness)
-
+ # Deal with any inf's or nan's
+ tf_at_f = mkinsertgap(pipeline, tf_at_f, replace_value = 0)
# Remove worthless data from computation of transfer function if we can
+ tf_at_f = mkgate(pipeline, tf_at_f, witness_downsampled[i], 1.0001 * witness_absent_value, attack_length = -((1.0 - filter_latency) * (filter_samples + resample_shift)), name = "linesub_witness_downsampled_gate_%d_%d_%d_%d" % (m, n, i, function_inst_num), queue_length = queue_length)
tf_at_f = mkgate(pipeline, tf_at_f, signal_exists, 1, attack_length = -((1.0 - filter_latency) * filter_samples), name = "linesub_signal_exists_gate_%d_%d_%d_%d" % (m, n, i, function_inst_num), queue_length = queue_length)
- tf_at_f = mkgate(pipeline, tf_at_f, witness_exists[i], 1, attack_length = -((1.0 - filter_latency) * (filter_samples + resample_shift)), name = "linesub_witness_exists_gate_%d_%d_%d_%d" % (m, n, i, function_inst_num), queue_length = queue_length)
if noisesub_gate_bit is not None:
tf_at_f = mkgate(pipeline, tf_at_f, noisesub_gate_bit, 1, attack_length = -((1.0 - filter_latency) * filter_samples), name = "linesub_gate_%d_%d_%d_%d" % (m, n, i, function_inst_num), queue_length = queue_length)
if any(amp_channels):
@@ -624,10 +628,14 @@ def remove_lines_with_witnesses(pipeline, signal, witnesses, freqs, freq_vars, f
tf_at_f = mkgate(pipeline, tf_at_f, amp_channels[m][n], 1e-35, attack_length = -((1.0 - filter_latency) * filter_samples), name = "linesub_amp_tf_gate_%d_%d_%d_%d" % (m, n, i, function_inst_num), queue_length = queue_length)
if(min_num_median_avg > 0 and (num_median > min_num_median_avg or num_avg > min_num_median_avg)):
tf_at_f = pipeparts.mktee(pipeline, tf_at_f)
- secular_change = detect_secular_change(pipeline, tf_at_f, (num_median + num_avg) / compute_rate, coherence_time = filter_length / 2, rate = compute_rate, secular_dur_min = (min(num_median, min_num_median_avg) + min(num_avg, min_num_median_avg)) / compute_rate, threshold = secular_change_threshold)
- secular_change = pipeparts.mkgeneric(pipeline, secular_change, "lal_property", update_samples = compute_rate * filter_length / 2, update_threshold = secular_change_threshold, name = "lal_property_secular_change_%d_%d_%d_%d" % (m, n, i, function_inst_num))
- tfs_at_f[i] = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_smoothkappas", no_default = True, reject_zeros = False, default_kappa_re = 0.0, min_array_size = min_num_median_avg, array_size = num_median, avg_array_size = num_avg, default_to_median = True, filter_latency = filter_latency)
- secular_change.connect("notify::current-average", reset_timestamped_secular_change, tfs_at_f[i], secular_change_threshold)
+ # Choose good parameters spaced approximately logarithmically
+ coherence_time = 2 * filter_length // 3
+ secular_dur = (num_median + num_avg) // compute_rate // 2
+ stochastic_dur = int(numpy.sqrt(secular_dur * coherence_time))
+ secular_change = detect_secular_change(pipeline, tf_at_f, secular_dur, stochastic_dur = stochastic_dur, coherence_time = coherence_time, rate = compute_rate)
+ secular_change = pipeparts.mkgeneric(pipeline, secular_change, "lal_property", update_samples = compute_rate * coherence_time, update_threshold = secular_change_threshold, update_after_gap = True, name = "lal_property_secular_change_%d_%d_%d_%d" % (m, n, i, function_inst_num))
+ tfs_at_f[i] = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_smoothkappas", no_default = True, reject_zeros = False, default_kappa_re = 0.0, min_array_size = min_num_median_avg, array_size = num_median, avg_array_size = num_avg, default_to_median = True, filter_latency = filter_latency, name = "linesub_smoothkappas_%d_%d_%d_%d" % (m, n, i, function_inst_num))
+ secular_change.connect("notify::timestamped-average", reset_secular_change, tfs_at_f[i])
else:
tfs_at_f[i] = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_smoothkappas", default_kappa_re = 0.0, array_size = num_median, avg_array_size = num_avg, default_to_median = True, filter_latency = filter_latency)
@@ -638,10 +646,11 @@ def remove_lines_with_witnesses(pipeline, signal, witnesses, freqs, freq_vars, f
if(i != j):
# Find transfer function between 2 witness channels at this frequency
tf_at_f = complex_division(pipeline, line_in_witnesses[j], line_in_witnesses[i])
-
+ # Deal with any inf's or nan's
+ tf_at_f = mkinsertgap(pipeline, tf_at_f, replace_value = 0)
# Remove worthless data from computation of transfer function if we can
- tf_at_f = mkgate(pipeline, tf_at_f, witness_exists[i], 1, attack_length = -((1.0 - filter_latency) * (filter_samples + resample_shift)), name = "linesub_witness_exists1_gate_%d_%d_%d_%d_%d" % (m, n, i, j, function_inst_num), queue_length = queue_length)
- tf_at_f = mkgate(pipeline, tf_at_f, witness_exists[j], 1, attack_length = -((1.0 - filter_latency) * (filter_samples + resample_shift)), name = "linesub_witness_exists2_gate_%d_%d_%d_%d_%d" % (m, n, i, j, function_inst_num), queue_length = queue_length)
+ tf_at_f = mkgate(pipeline, tf_at_f, witness_downsampled[i], 1.0001 * witness_absent_value, attack_length = -((1.0 - filter_latency) * (filter_samples + resample_shift)), name = "linesub_witness_downsampled1_gate_%d_%d_%d_%d_%d" % (m, n, i, j, function_inst_num), queue_length = queue_length)
+ tf_at_f = mkgate(pipeline, tf_at_f, witness_downsampled[j], 1.0001 * witness_absent_value, attack_length = -((1.0 - filter_latency) * (filter_samples + resample_shift)), name = "linesub_witness_downsampled2_gate_%d_%d_%d_%d_%d" % (m, n, i, j, function_inst_num), queue_length = queue_length)
if noisesub_gate_bit is not None:
tf_at_f = mkgate(pipeline, tf_at_f, noisesub_gate_bit, 1, attack_length = -((1.0 - filter_latency) * filter_samples), name = "linesub_gate_%d_%d_%d_%d_%d" % (m, n, i, j, function_inst_num), queue_length = queue_length)
if any(amp_channels):
@@ -649,10 +658,14 @@ def remove_lines_with_witnesses(pipeline, signal, witnesses, freqs, freq_vars, f
tf_at_f = mkgate(pipeline, tf_at_f, amp_channels[m][n], 1e-35, attack_length = -((1.0 - filter_latency) * filter_samples), name = "linesub_amp_tf_gate_%d_%d_%d_%d_%d" % (m, n, i, j, function_inst_num), queue_length = queue_length)
if(min_num_median_avg > 0 and (num_median > min_num_median_avg or num_avg > min_num_median_avg)):
tf_at_f = pipeparts.mktee(pipeline, tf_at_f)
- secular_change = detect_secular_change(pipeline, tf_at_f, (num_median + num_avg) / compute_rate, coherence_time = filter_length / 2, rate = compute_rate, secular_dur_min = (min(num_median, min_num_median_avg) + min(num_avg, min_num_median_avg)) / compute_rate, threshold = secular_change_threshold)
- secular_change = pipeparts.mkgeneric(pipeline, secular_change, "lal_property", update_samples = compute_rate * filter_length / 2, update_threshold = secular_change_threshold, name = "lal_property_secular_change_%d_%d_%d_%d_%d" % (m, n, i, j, function_inst_num))
+ # Choose good parameters spaced approximately logarithmically
+ coherence_time = 2 * filter_length // 3
+ secular_dur = (num_median + num_avg) // compute_rate // 2
+ stochastic_dur = int(numpy.sqrt(secular_dur * coherence_time))
+ secular_change = detect_secular_change(pipeline, tf_at_f, secular_dur, stochastic_dur = stochastic_dur, coherence_time = coherence_time, rate = compute_rate)
+ secular_change = pipeparts.mkgeneric(pipeline, secular_change, "lal_property", update_samples = compute_rate * coherence_time, update_threshold = secular_change_threshold, update_after_gap = True, name = "lal_property_secular_change_%d_%d_%d_%d_%d" % (m, n, i, j, function_inst_num))
tfs_at_f[(i + 1) * len(witnesses[m]) + j] = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_smoothkappas", no_default = True, reject_zeros = False, default_kappa_re = 0.0, min_array_size = min_num_median_avg, array_size = num_median, avg_array_size = num_avg, default_to_median = True, filter_latency = filter_latency)
- secular_change.connect("notify::current-average", reset_timestamped_secular_change, tfs_at_f[(i + 1) * len(witnesses[m]) + j], secular_change_threshold)
+ secular_change.connect("notify::timestamped-average", reset_secular_change, tfs_at_f[(i + 1) * len(witnesses[m]) + j])
else:
tfs_at_f[(i + 1) * len(witnesses[m]) + j] = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_smoothkappas", default_kappa_re = 0.0, array_size = num_median, avg_array_size = num_avg, default_to_median = True, filter_latency = filter_latency)
@@ -1896,25 +1909,9 @@ def update_timestamped_property(prop_maker, arg, prop_taker, maker_prop_name, ta
prop = prop_maker.get_property(maker_prop_name)
cs.set(int(prop[0] * Gst.SECOND), prefactor * prop[1])
-def reset_secular_change(prop_maker, arg, prop_taker, threshold):
- prop = prop_maker.get_property("current_average")
- if float(prop) > threshold:
- prop_taker.set_property("reset", True)
-
-def reset_timestamped_secular_change(prop_maker, arg, prop_taker, threshold):
-
- binding = prop_taker.get_control_binding("reset")
- if(binding is not None):
- cs = binding.get_property('control_source')
- else:
- cs = GstController.InterpolationControlSource.new()
- cs.set_property('mode', GstController.InterpolationMode.NONE) # no interpolation
- binding = GstController.DirectControlBinding.new_absolute(prop_taker, "reset", cs)
- prop_taker.add_control_binding(binding)
-
+def reset_secular_change(prop_maker, arg, prop_taker):
prop = prop_maker.get_property("timestamped_average")
- if float(prop[1]) > threshold:
- cs.set(int(prop[0] * Gst.SECOND), True)
+ prop_taker.set_property("reset_time", int(2 + prop[0]) * Gst.SECOND)
def update_SS_demod_freqs(prop_maker, arg, prop_taker, maker_prop_name, taker_prop_name, prefactor):
@@ -2615,7 +2612,7 @@ def detect_secular_change(pipeline, head, secular_dur, stochastic_dur = 0, coher
nongap = pipeparts.mkcapsfilter(pipeline, nongap, "audio/x-raw,format=U32LE,rate=%d,channels=1,channel-mask=(bitmask)0x0" % rate)
nongap = pipeparts.mkgeneric(pipeline, nongap, "lal_typecast")
nongap = pipeparts.mkcapsfilter(pipeline, nongap, "audio/x-raw,format=F64LE,rate=%d,channels=1,channel-mask=(bitmask)0x0" % rate)
- head = pipeparts.mktee(pipeline, mkinsertgap(pipeline, head, insert_gap = False, remove_gap = True))
+ head = pipeparts.mktee(pipeline, head)
sqrt_num = pipeparts.mkgeneric(pipeline, head, "lal_stdev", mode = 2, array_size = int(secular_dur / coherence_time), coherence_length = int(coherence_time * rate), min_array_size = int(secular_dur_min / coherence_time), percentile = percentile)
num = mkpow(pipeline, sqrt_num, exponent = 2.0)
den = pipeparts.mkgeneric(pipeline, head, "lal_stdev", mode = 2, array_size = int(stochastic_dur / coherence_time), coherence_length = int(coherence_time * rate), percentile = percentile)
@@ -2623,14 +2620,17 @@ def detect_secular_change(pipeline, head, secular_dur, stochastic_dur = 0, coher
den = pipeparts.mkgeneric(pipeline, den, "lal_smoothkappas", array_size = secular_dur * rate, avg_array_size = 1, default_kappa_re = 0, default_to_median = True, no_default = True, min_array_size = secular_dur_min * rate, reject_zeros = False)
inv_den = complex_inverse(pipeline, den)
ratio = mkmultiplier(pipeline, [num, inv_den, nongap])
- ratio = mkinsertgap(pipeline, ratio, insert_gap = False, replace_value = 0)
+ # Get rid of inf's and nan's
+ ratio = mkinsertgap(pipeline, ratio, replace_value = 1e49, insert_gap = False)
+ # Add GAP flag back if it was removed
+ ratio = mkinsertgap(pipeline, ratio, bad_data_intervals = [-1e50, 0, 0, 1e50], replace_value = 0)
if secular_dur_min > 0 and secular_dur_min < secular_dur and threshold > 1:
# Reset the calculation anytime ratio > threshold, which indicates secular change
ratio = pipeparts.mktee(pipeline, ratio)
- ratio_prop = pipeparts.mkgeneric(pipeline, ratio, "lal_property", update_samples = coherence_time * rate, update_threshold = threshold)
- ratio_prop.connect("notify::current-average", reset_timestamped_secular_change, sqrt_num, threshold)
- ratio_prop.connect("notify::current-average", reset_timestamped_secular_change, den, threshold)
+ ratio_prop = pipeparts.mkgeneric(pipeline, ratio, "lal_property", update_samples = coherence_time * rate, update_threshold = threshold, update_after_gap = True)
+ ratio_prop.connect("notify::timestamped-average", reset_secular_change, sqrt_num, secular_dur_min)
+ ratio_prop.connect("notify::timestamped-average", reset_secular_change, den, secular_dur_min)
return ratio
diff --git a/tests/check_calibration/Makefile b/tests/check_calibration/Makefile
index 2ad4c28cd5b42f905f22a570cf1c3ebb594522fb..1a58247c14d0a2b2935c617932f0fc394eeeaf57 100644
--- a/tests/check_calibration/Makefile
+++ b/tests/check_calibration/Makefile
@@ -4,13 +4,14 @@
#################################
# which interferometer (H or L)
-IFO = L
+IFO = H
# determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14)
OBSRUN = O3
# Determines whether to find calibration filters in the old SVN or the git repo. Starting with PreER15 and O4, we switched to the git repo. Set to 'svn' or 'git'.
FILTERSRC = svn
-START = $(shell echo 1406743118 | bc)
+START = $(shell echo 1409576018 | bc)
+# 1409576018 or 1410782418 gstlal-calibration-1.5.8 linesub issue
# 1399057926 line sub not working
# 1371862818 oversubtraction study
# 1266880202 Comb
@@ -20,7 +21,8 @@ START = $(shell echo 1406743118 | bc)
# 1238288418
# 1269090018
# 1269114645-21449
-END = $(shell echo 1406748618 | bc)
+END = $(shell echo 1409604100 | bc)
+# 1409604100 or 1410839118 gstlal-calibration-1.5.8 linesub issue
# 1399086062 line sub not working
# 1371911734 oversubtraction study
# 1266888034 Comb
@@ -421,7 +423,7 @@ noise_subtraction_ASD_GDS: $(IFO)1_hoft_GDS_frames.cache $(IFO)1_hoft_GDS_frames
python3 plot_ASD.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache-list $(IFO)1_hoft_GDS_frames.cache,$(IFO)1_hoft_GDS_frames_shortTF.cache,$(IFO)1_hoft_GDS_frames_longTF.cache,$(IFO)1_hoft_GDS_frames.cache --channel-list GDS-CALIB_STRAIN,GDS-CALIB_STRAIN_NOLINES,GDS-CALIB_STRAIN_NOLINES,GDS-CALIB_STRAIN_NOLINES --sample-rate 16384 --fft-time 512 --fft-spacing 256 --window=0 --freq-res 0.05 --frequency-min 1 --frequency-max=8192 --ASD-max=1e-17 --labels 'strain,128s\ estimate,4096s\ estimate,switch'
noise_subtraction_spectrogram_GDS: $(IFO)1_hoft_GDS_frames.cache
- python3 plot_spectrogram.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache $(IFO)1_hoft_GDS_frames.cache --channel-name GDS-CALIB_STRAIN_NOLINES --sample-rate 16384 --fft-time 1 --fft-spacing 0.25 --window=3 --frequency-min 10 --frequency-max=1000 --update-time 4 --amplitude-min 0.01 --amplitude-max 100
+ python3 plot_spectrogram.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache $(IFO)1_hoft_GDS_frames.cache --channel-name GDS-CALIB_STRAIN_NOLINES --sample-rate 16384 --window=3 --frequency-min 10 --frequency-max=1000 --amplitude-min 0.01 --amplitude-max 100
noise_subtraction_ASD_GDS1: $(IFO)1_hoft_GDS_frames.cache $(IFO)1_hoft_GDS_frames_shortTF.cache $(IFO)1_hoft_GDS_frames_longTF.cache
python3 plot_ASD.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache-list $(IFO)1_hoft_GDS_frames.cache,$(IFO)1_hoft_GDS_frames_shortTF.cache,$(IFO)1_hoft_GDS_frames_longTF.cache,$(IFO)1_hoft_GDS_frames.cache --channel-list GDS-CALIB_STRAIN,GDS-CALIB_STRAIN_NOLINES,GDS-CALIB_STRAIN_NOLINES,GDS-CALIB_STRAIN_NOLINES --sample-rate 16384 --fft-time 4096 --fft-spacing 2048 --freq-res 0.001 --frequency-min 410.27 --frequency-max=410.33 --ASD-max=1e-19 --labels 'strain,128s\ estimate,4096s\ estimate,switch'
diff --git a/tests/check_calibration/plot_spectrogram.py b/tests/check_calibration/plot_spectrogram.py
index 76c64eea6340b043b108d1249d43fb5ab80113dd..3f084b9b75b292d676c9be0518422485f05a4f7b 100644
--- a/tests/check_calibration/plot_spectrogram.py
+++ b/tests/check_calibration/plot_spectrogram.py
@@ -67,11 +67,11 @@ parser.add_option("--ifo", metavar = "name", type = str, help = "Name of the int
parser.add_option("--frame-cache", metavar = "name", type = str, help = "Frame cache file that contains data from which to compute the spectrogram")
parser.add_option("--channel-name", metavar = "name", type = str, default = None, help = "Names of channel that contains data from which to compute the spectrogram")
parser.add_option("--sample-rate", metavar = "Hz", type = str, default = '16384', help = "Sample rate of input data.")
-parser.add_option("--fft-time", metavar = "seconds", type = float, default = 1, help = "Duration in seconds of FFTs used to compute spectrogram")
-parser.add_option("--fft-spacing", metavar = "seconds", type = float, default = 0.25, help = "Spacing in seconds of FFTs used to compute spectrogram")
+parser.add_option("--fft-time", metavar = "seconds", type = float, default = 4, help = "Duration in seconds of FFTs used to compute spectrogram")
+parser.add_option("--fft-spacing", metavar = "seconds", type = float, default = None, help = "Spacing in seconds of FFTs used to compute spectrogram. Default is auto.")
parser.add_option("--window", type = int, default = 3, help = "Which window function to use to window the data when taking FFTs. Options are 0 (DPSS), 1 (kaiser), 2 (Dolph-Chebyshev), 3 (Blackman, default), 4 (Hann), 5 (None)")
parser.add_option("--freq-res", metavar = "Hz", type = float, default = 1.0, help = "Frequency resolution of spectrogram. Only applies when using Kaiser, DPSS, and Dolph Chebyshev windows.")
-parser.add_option("--update-time", type = float, default = 4, help = "Amount of time between consecutive ASDs on the spectrogram")
+parser.add_option("--update-time", type = float, default = None, help = "Amount of time between consecutive ASDs on the spectrogram. Default is auto.")
parser.add_option("--frequency-min", type = float, default = 10, help = "Minimum frequency for plot")
parser.add_option("--frequency-max", type = float, default = 8192, help = "Maximum frequency for plot")
parser.add_option("--ASD-min", type = float, default = 1e-24, help = "Minimum for ASD plot")
@@ -88,6 +88,21 @@ frame_cache = options.frame_cache
channel = options.channel_name
chan_list = [(ifo, channel)]
sr = sample_rate = int(options.sample_rate)
+update_time = options.update_time
+dur = options.gps_end_time - options.gps_start_time
+if update_time is None:
+ # ~1000 samples is about right
+ update_time = dur / 1000.0
+ # Make it a power of 2
+ update_time = 2**int(round(np.log2(update_time)))
+ # It should be more than the fft_time
+ while update_time < options.fft_time:
+ update_time *= 2
+fft_spacing = options.fft_spacing
+if fft_spacing is None:
+ fft_spacing = update_time / 8
+ while fft_spacing < options.fft_time / 4 and fft_spacing < update_time:
+ fft_spacing *= 2
#
# =============================================================================
@@ -112,9 +127,9 @@ def compute_spectrogram(pipeline, name):
data = calibration_parts.mkresample(pipeline, data, 5, False, sr)
data = pipeparts.mktee(pipeline, data)
# Compute ASDs for the spectrogram
- pipeparts.mkgeneric(pipeline, data, "lal_asd", fft_samples = int(options.fft_time * sr), overlap_samples = int((options.fft_time - options.fft_spacing) * sr), window_type = int(options.window), frequency_resolution = float(options.freq_res), update_time = int(options.update_time), use_td_median = True, filename = '%s_%s_spectrogram_%d-%d.txt' % (ifo, channel, options.gps_start_time, options.gps_end_time - options.gps_start_time))
+ pipeparts.mkgeneric(pipeline, data, "lal_asd", fft_samples = int(options.fft_time * sr), overlap_samples = int((options.fft_time - fft_spacing) * sr), window_type = int(options.window), frequency_resolution = float(options.freq_res), update_time = update_time, use_td_median = True, filename = '%s_%s_spectrogram_%d-%d.txt' % (ifo, channel, options.gps_start_time, dur))
# Compute one median ASD to normalize the above ASDs
- pipeparts.mkgeneric(pipeline, data, "lal_asd", fft_samples = int(options.fft_time * sr), overlap_samples = int((options.fft_time - options.fft_spacing) * sr), window_type = int(options.window), frequency_resolution = float(options.freq_res), update_time = int(options.gps_end_time - options.gps_start_time), use_td_median = True, filename = '%s_%s_ASD_%d-%d.txt' % (ifo, channel, options.gps_start_time, options.gps_end_time - options.gps_start_time))
+ pipeparts.mkgeneric(pipeline, data, "lal_asd", fft_samples = int(options.fft_time * sr), overlap_samples = int((options.fft_time - fft_spacing) * sr), window_type = int(options.window), frequency_resolution = float(options.freq_res), update_time = 2 * int(dur), use_td_median = True, filename = '%s_%s_ASD_%d-%d.txt' % (ifo, channel, options.gps_start_time, dur))
#
# done
@@ -135,20 +150,20 @@ def compute_spectrogram(pipeline, name):
test_common.build_and_run(compute_spectrogram, "compute_spectrogram", segment = segments.segment((LIGOTimeGPS(0, 1000000000 * options.gps_start_time), LIGOTimeGPS(0, 1000000000 * options.gps_end_time))))
# Unnormalized spectrogram data
-spec = np.loadtxt('%s_%s_spectrogram_%d-%d.txt' % (ifo, channel, options.gps_start_time, options.gps_end_time - options.gps_start_time))
+spec = np.loadtxt('%s_%s_spectrogram_%d-%d.txt' % (ifo, channel, options.gps_start_time, dur))
# ASD to normalize it
-asd = np.loadtxt('%s_%s_ASD_%d-%d.txt' % (ifo, channel, options.gps_start_time, options.gps_end_time - options.gps_start_time))
+asd = np.loadtxt('%s_%s_ASD_%d-%d.txt' % (ifo, channel, options.gps_start_time, dur))
# Frequency vector
fvec = np.arange(0, sr / 2.0 + sr / 4.0 / len(asd), sr / 2.0 / (len(asd) - 1))
# Time vector
num_t = len(spec) // len(asd)
-max_t = options.update_time * (num_t - 1)
-tvec = np.arange(0, max_t + options.update_time / 2, options.update_time)
+max_t = update_time * (num_t - 1)
+tvec = np.arange(0, max_t + update_time / 2, update_time)
-t0 = options.gps_start_time + options.update_time / 2
+t0 = options.gps_start_time + update_time / 2
# Normalize spectrogram
asd = np.transpose(asd)[1]
@@ -189,7 +204,7 @@ plt.title(r'${\rm %s:%s \ Spectrogram}$' % (ifo, channel.replace('_', '\_')))
ticks_and_grid(ax, ymin = options.frequency_min, ymax = options.frequency_max, yscale = freq_scale, xscale = 'linear')
-plt.savefig('%s_%s_spectrogram_%d-%dHz_%d-%d.png' % (ifo, channel, int(options.frequency_min), int(options.frequency_max), options.gps_start_time, options.gps_end_time - options.gps_start_time))
+plt.savefig('%s_%s_spectrogram_%d-%dHz_%d-%d.png' % (ifo, channel, int(options.frequency_min), int(options.frequency_max), options.gps_start_time, dur))
plt.close()
# Median ASD for reference
@@ -201,6 +216,6 @@ plt.xlabel(r'${\rm Frequency \ [Hz]}$')
ticks_and_grid(plt.gca(), xmin = options.frequency_min, xmax = options.frequency_max, ymin = options.ASD_min, ymax = options.ASD_max, xscale = freq_scale, yscale = ASD_scale)
plt.tight_layout()
-plt.savefig('%s_%s_ASD_%d-%dHz_%d-%d.png' % (ifo, channel, int(options.frequency_min), int(options.frequency_max), options.gps_start_time, options.gps_end_time - options.gps_start_time))
+plt.savefig('%s_%s_ASD_%d-%dHz_%d-%d.png' % (ifo, channel, int(options.frequency_min), int(options.frequency_max), options.gps_start_time, dur))
diff --git a/tests/tests_pytest/STV_error.py b/tests/tests_pytest/STV_error.py
index 9c2374cbb6ebd99f47387f39bfbb065e9d9fbcac..374f58d2d9cdaed3932238ab23fb9d19ee0f0d02 100644
--- a/tests/tests_pytest/STV_error.py
+++ b/tests/tests_pytest/STV_error.py
@@ -20,7 +20,7 @@ import numpy as np
def check_bits(filename):
# Numbers corresponding to each bit
bits = []
- for i in range(21):
+ for i in range(20):
bits.append(2**i)
# Bit definitions
@@ -45,7 +45,7 @@ def check_bits(filename):
bitnames.append("line_sub_bit")
bitnames.append("noise_sub_bit")
bitnames.append("noise_sub_gate_bit")
- bitnames.append("nonsens_sub_bit")
+ #bitnames.append("nonsens_sub_bit")
standard = np.loadtxt('tests/tests_pytest/State_Vector_data/%s_standard.txt' % filename)
data = np.loadtxt('tests/tests_pytest/State_Vector_data/%s.txt' % filename)
@@ -96,7 +96,7 @@ def check_bits(filename):
assert consecutive_failures < 160, "%s failure!" % bitnames[i] # 10 seconds at a 16-Hz sample rate
e_file.write("%s passed\n" % bitnames[i])
e_file.close()
- for i in [17, 18, 20]:
+ for i in [17, 18]:
# Require agreement at a majority of samples
e_file = open('tests/tests_pytest/error_results.txt', 'a')
e_file.write("Testing %s %s\n" % (filename, bitnames[i]))
diff --git a/tests/tests_pytest/error.py b/tests/tests_pytest/error.py
index 1c46e6698bf8f349eda2cab3f188b8ee3104d71c..4eeb447be708f64ceddc8e6d5ef1b26b6dd474e4 100644
--- a/tests/tests_pytest/error.py
+++ b/tests/tests_pytest/error.py
@@ -143,5 +143,8 @@ def rms(d_type, typ='A'):
e_file.close()
print(s2, file = sys.stderr)
- assert rms < 0.002, "%s: error too large!" % name[i]
+ if 'clean hoft' in name[i]:
+ assert rms < 0.01, "%s: error too large!" % name[i]
+ else:
+ assert rms < 0.002, "%s: error too large!" % name[i]
diff --git a/tests/tests_pytest/filters/gstlal_compute_strain_ApproxLaterStart_H1.ini b/tests/tests_pytest/filters/gstlal_compute_strain_ApproxLaterStart_H1.ini
index 0c7a2a42753a9dd97ca1f4e7fe0eecb973c5e331..45801a472773827848b6292314b0a9f5e96c89be 100644
--- a/tests/tests_pytest/filters/gstlal_compute_strain_ApproxLaterStart_H1.ini
+++ b/tests/tests_pytest/filters/gstlal_compute_strain_ApproxLaterStart_H1.ini
@@ -398,9 +398,11 @@ CBCHWInjOffBitmask: 9
BurstHWInjOffBitmask: 17
DetCharHWInjOffBitmask: 33
StochHWInjOffBitmask: 65
-NoiseSubGateBitmask: 512
+#NoiseSubGateBitmask: 512
# If we require NoiseSubGateChannel > some value, use this instead. Otherwise, set to "None"
-NoiseSubGateThreshold: 599
+#NoiseSubGateThreshold: 599
+# Require the noise subtraction gate channel to be equal to this value for the gate to be open
+NoiseSubGateValue: 600
FilterClockBitmaskList: 1
UndisturbedOKBitmask: 1
diff --git a/tests/tests_pytest/filters/gstlal_compute_strain_Approx_H1.ini b/tests/tests_pytest/filters/gstlal_compute_strain_Approx_H1.ini
index 7fd0e7c8b3d2a5941b784ef65abe940b13b05ab3..345b9ca46e885027bf43d326ae51ff00c6d73575 100644
--- a/tests/tests_pytest/filters/gstlal_compute_strain_Approx_H1.ini
+++ b/tests/tests_pytest/filters/gstlal_compute_strain_Approx_H1.ini
@@ -398,9 +398,11 @@ CBCHWInjOffBitmask: 9
BurstHWInjOffBitmask: 17
DetCharHWInjOffBitmask: 33
StochHWInjOffBitmask: 65
-NoiseSubGateBitmask: 512
+#NoiseSubGateBitmask: 512
# If we require NoiseSubGateChannel > some value, use this instead. Otherwise, set to "None"
-NoiseSubGateThreshold: 599
+#NoiseSubGateThreshold: 599
+# Require the noise subtraction gate channel to be equal to this value for the gate to be open
+NoiseSubGateValue: 600
FilterClockBitmaskList: 1
UndisturbedOKBitmask: 1
diff --git a/tests/tests_pytest/filters/gstlal_compute_strain_Exact_H1.ini b/tests/tests_pytest/filters/gstlal_compute_strain_Exact_H1.ini
index e62770745587cf7c42f7a7c5b73e561a52a52fbd..9de4370a20591a80332f0c6ecae62392dcd76345 100644
--- a/tests/tests_pytest/filters/gstlal_compute_strain_Exact_H1.ini
+++ b/tests/tests_pytest/filters/gstlal_compute_strain_Exact_H1.ini
@@ -398,9 +398,11 @@ CBCHWInjOffBitmask: 9
BurstHWInjOffBitmask: 17
DetCharHWInjOffBitmask: 33
StochHWInjOffBitmask: 65
-NoiseSubGateBitmask: 512
+#NoiseSubGateBitmask: 512
# If we require NoiseSubGateChannel > some value, use this instead. Otherwise, set to "None"
-NoiseSubGateThreshold: 599
+#NoiseSubGateThreshold: 599
+# Require the noise subtraction gate channel to be equal to this value for the gate to be open
+NoiseSubGateValue: 600
FilterClockBitmaskList: 1
UndisturbedOKBitmask: 1
diff --git a/tests/tests_pytest/latency.py b/tests/tests_pytest/latency.py
index 1495f721cbdda40c0d4fdb53eb5a8054e1615258..423cec3b818a2189b7ee3a92e23d1ab6bc439fa5 100644
--- a/tests/tests_pytest/latency.py
+++ b/tests/tests_pytest/latency.py
@@ -35,26 +35,31 @@ def Latency():
i_lst = [list(i_gps).index(x) for x in list(np.intersect1d(i_gps,o_gps))] #sub a and b for left and right collumn of output, returns a list of the indexs of a to keep
lat = o_real - i_real[i_lst]
- # Allow 100 s for pipeline startup processes to settle, and remove the last frame, which may be a partial frame
- settled_lat = lat[100:-1]
+ # Allow 200 s for pipeline startup processes to settle, and remove the last frame, which may be a partial frame
+ settled_lat = lat[50:-1]
# Check for an increasing trend
var = settled_lat - np.mean(settled_lat)
+ # Take a running 20th percentile to reduce the impact of outliers
+ percentile_length = 90 # 6 minutes
+ percentile_var = np.zeros(len(var) - percentile_length + 1)
+ for i in range(len(percentile_var)):
+ percentile_var[i] = np.percentile(var[i : i + percentile_length], 20)
# lat_increase estimates the latency increase in seconds from start to finish
lat_increase = 0.0
- for i in range(len(var)):
- lat_increase += var[i] * 12 * (i - (len(var) - 1) / 2.0) / len(var) / len(var)
+ for i in range(len(percentile_var)):
+ lat_increase += percentile_var[i] * 12 * (i - (len(percentile_var) - 1) / 2.0) / len(percentile_var) / len(percentile_var)
e_file = open('tests/tests_pytest/error_results.txt', 'a')
e_file.write("Mean latency = %fs\n" % np.mean(settled_lat))
e_file.write("Max latency = %fs\n" % max(settled_lat))
e_file.write("%d of %d frames with latency over 6s\n" % (sum(i > 6 for i in settled_lat), len(settled_lat)))
- e_file.write("Latency increased by %fs over %fs of data\n" % (lat_increase, o_gps[-1] - o_gps[0]))
+ e_file.write("Latency increased by %fs over %fs of data\n" % (lat_increase, 4 * (len(percentile_var) - 1)))
e_file.close()
print("Mean latency = {}s".format(np.mean(settled_lat)), file = sys.stderr)
print("Max latency = {}s".format(max(settled_lat)), file = sys.stderr)
print("{} of {} frames with latency over 6s".format(sum(i > 6 for i in settled_lat), len(settled_lat)), file = sys.stderr)
- print("Latency increased by {}s over {}s of data".format(lat_increase, o_gps[-1] - o_gps[0]), file = sys.stderr)
+ print("Latency increased by {}s over {}s of data".format(lat_increase, 4 * (len(percentile_var) - 1)), file = sys.stderr)
np.savetxt("tests/tests_pytest/latency.txt", np.transpose(np.array([o_gps, lat])))
plot_latency(o_gps, lat)
@@ -62,13 +67,13 @@ def Latency():
# 4-s latency is typical in this configuration. More than 5 s could
# indicate a problem, while less than 3 s could mean that the CI
# latency-testing pipeline is not running as it should be.
- assert np.mean(settled_lat) < 5
- assert np.mean(settled_lat) > 3
+ assert np.mean(settled_lat) < 5, "Mean latency too large (%f s > 5 s)" % np.mean(settled_lat)
+ assert np.mean(settled_lat) > 3, "Mean latency too small (%f s < 3 s)" % np.mean(settled_lat)
# There should be no large excursions (This may occasionally fail)
- assert max(settled_lat) < 10
- assert min(settled_lat) > 0
+ assert max(settled_lat) < 10, "Max latency too large (%f s > 10 s)" % max(settled_lat)
+ assert min(settled_lat) > 0, "Min latency too small (%f s < 0 s)" % min(settled_lat)
# There should not be numerous small excursions
- assert sum(i > 6 for i in settled_lat) < 5
- assert sum(i < 2 for i in settled_lat) < 5
+ assert sum(i > 6 for i in settled_lat) < 10, "Too many latency excursions (%d frames had latency > 6 s)" % sum(i > 6 for i in settled_lat)
+ assert sum(i < 2 for i in settled_lat) < 10, "Too many latency excursions (%d frames had latency < 2 s)" % sum(i < 2 for i in settled_lat)
assert lat_increase < 0.5
diff --git a/tests/tests_pytest/run_calib_pipeline.py b/tests/tests_pytest/run_calib_pipeline.py
index 70d21239c3e0242f6a23d3eb73efd6e6c4e31796..be3ccb66cc2efe15ba0fc9b7983ebcecc6fc12d9 100644
--- a/tests/tests_pytest/run_calib_pipeline.py
+++ b/tests/tests_pytest/run_calib_pipeline.py
@@ -15,6 +15,7 @@
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
import os
+from threading import Thread
gps_start_time = 1404304216
@@ -30,38 +31,38 @@ later_frame_start = later_start_time + 32 + (4 - later_start_time % 4) % 4
expected_fnum = (gps_end_time - gps_frame_start) // 2 + (gps_end_time - later_frame_start) // 4
-def Run_calib():
-
- os.system("ls tests/tests_pytest/frames/raw/H-H1_R-*.gwf | lal_path2cache > tests/tests_pytest/raw_frames.cache")
- os.system("mkdir -p tests/tests_pytest/frames/GDS")
-
- e_file = open('tests/tests_pytest/error_results.txt', 'a')
- e_file.write('Running Approx-TDCFs calibration pipeline\n')
- e_file.close()
+def Run_calib_approx_exact():
os.system("GST_DEBUG=3 gstlal_compute_strain --gps-start-time %d --gps-end-time %d --frame-cache tests/tests_pytest/raw_frames.cache --output-path tests/tests_pytest/frames/GDS --frame-duration=4 --frames-per-file=1 --wings=0 --config-file tests/tests_pytest/filters/gstlal_compute_strain_Approx_H1.ini --filters-file-name tests/tests_pytest/filters/gstlal_compute_strain_C00_filters_H1_20240330T211519Z.npz" % (gps_start_time, gps_end_time))
- e_file = open('tests/tests_pytest/error_results.txt', 'a')
- e_file.write('Finished running Approx-TDCFs calibration pipeline\n')
- e_file.close()
os.system("ls tests/tests_pytest/frames/GDS/H-H1GDS_Approx-*.gwf | lal_path2cache > tests/tests_pytest/GDS_Approx_frames.cache")
- e_file = open('tests/tests_pytest/error_results.txt', 'a')
- e_file.write('Running Exact-TDCFs calibration pipeline\n')
- e_file.close()
os.system("GST_DEBUG=3 gstlal_compute_strain --gps-start-time %d --gps-end-time %d --frame-cache tests/tests_pytest/raw_frames.cache --output-path tests/tests_pytest/frames/GDS --frame-duration=4 --frames-per-file=1 --wings=0 --config-file tests/tests_pytest/filters/gstlal_compute_strain_Exact_H1.ini --filters-file-name tests/tests_pytest/filters/gstlal_compute_strain_C00_filters_H1_20240330T211519Z.npz" % (gps_start_time, gps_end_time))
- e_file = open('tests/tests_pytest/error_results.txt', 'a')
- e_file.write('Finished running Exact-TDCFs calibration pipeline\n')
- e_file.close()
os.system("ls tests/tests_pytest/frames/GDS/H-H1GDS_Exact-*.gwf | lal_path2cache > tests/tests_pytest/GDS_Exact_frames.cache")
- e_file = open('tests/tests_pytest/error_results.txt', 'a')
- e_file.write('Running Approx-TDCFs calibration pipeline with a later start time\n')
- e_file.close()
+
+def Run_calib_later_start():
os.system("GST_DEBUG=3 gstlal_compute_strain --gps-start-time %d --gps-end-time %d --frame-cache tests/tests_pytest/raw_frames.cache --output-path tests/tests_pytest/frames/GDS --frame-duration=4 --frames-per-file=1 --wings=0 --config-file tests/tests_pytest/filters/gstlal_compute_strain_ApproxLaterStart_H1.ini --filters-file-name tests/tests_pytest/filters/gstlal_compute_strain_C00_filters_H1_20240330T211519Z.npz" % (later_start_time, gps_end_time))
+ os.system("ls tests/tests_pytest/frames/GDS/H-H1GDS_LaterStart-*.gwf | lal_path2cache > tests/tests_pytest/GDS_LaterStart_frames.cache")
+
+
+def Run_calib():
+ os.system("ls tests/tests_pytest/frames/raw/H-H1_R-*.gwf | lal_path2cache > tests/tests_pytest/raw_frames.cache")
+ os.system("mkdir -p tests/tests_pytest/frames/GDS")
+
e_file = open('tests/tests_pytest/error_results.txt', 'a')
- e_file.write('Finished running Approx-TDCFs calibration pipeline with a later start time\n')
+ e_file.write('Running calibration pipelines\n')
e_file.close()
- os.system("ls tests/tests_pytest/frames/GDS/H-H1GDS_LaterStart-*.gwf | lal_path2cache > tests/tests_pytest/GDS_LaterStart_frames.cache")
+
+ approx_exact_thread = Thread(target = Run_calib_approx_exact)
+ later_start_thread = Thread(target = Run_calib_later_start)
+ approx_exact_thread.start()
+ later_start_thread.start()
+ approx_exact_thread.join()
+ later_start_thread.join()
fnum = len(os.listdir('tests/tests_pytest/frames/GDS'))
- assert fnum == expected_fnum
+ e_file = open('tests/tests_pytest/error_results.txt', 'a')
+ e_file.write('Finished running calibration pipelines\n')
+ e_file.write('%d frames expected from calibration pipelines\n%d frames produced by calibration pipelines\n' % (expected_fnum, fnum))
+ e_file.close()
+ assert fnum == expected_fnum, "Calibration pipelines should produce %d frames but only produced %d." % (expected_fnum, fnum)
diff --git a/tests/tests_pytest/update_standards.py b/tests/tests_pytest/update_standards.py
new file mode 100644
index 0000000000000000000000000000000000000000..f5a5d3c5b05a0752ed95ddc79d89716dd2a675a4
--- /dev/null
+++ b/tests/tests_pytest/update_standards.py
@@ -0,0 +1,111 @@
+
+import numpy as np
+import os
+
+
+update_act_standard_files()
+update_pcal_standard_files()
+update_ASD_standard_files()
+update_TDCF_standard_files()
+update_STV_standard_files()
+
+
+def update_act_standard_files():
+ act_data_files = os.popen("ls act_data/*exc_0*.txt").read().split('\n')[:-1] + os.popen("ls act_data/*exc_1*.txt").read().split('\n')[:-1]
+ for act_data_file in act_data_files:
+ act_data_standard_file = act_data_file.split('exc')[0] + 'exc_standard' + act_data_file.split('exc')[1]
+ act_data_standard = np.loadtxt(act_data_standard_file)
+ act_data_standard_times = np.transpose(act_data_standard)[0]
+ act_data = np.loadtxt(act_data_file)
+ act_data_standard = np.zeros([len(act_data_standard), 3])
+ i = 0
+ successes = 0
+ for j in range(len(act_data_standard_times)):
+ while act_data_standard_times[j] != act_data[i][0] and i < len(act_data):
+ i += 1
+ if act_data_standard_times[j] == act_data[i][0]:
+ act_data_standard[j] = act_data[i]
+ successes += 1
+ if successes != len(act_data_standard_times):
+ print("%s: mismatch between data in new file and data in standard file. Expected %d timestamp matches, found %d" % (act_data_file, len(act_data_standard_times), successes))
+ else:
+ os.system("rm %s" % act_data_standard_file)
+ np.savetxt(act_data_standard_file, act_data_standard)
+
+
+def update_pcal_standard_files():
+ pcal_data_files = os.popen("ls pcal_data/*Pcal_0*.txt").read().split('\n')[:-1] + os.popen("ls pcal_data/*Pcal_1*.txt").read().split('\n')[:-1]
+ for pcal_data_file in pcal_data_files:
+ pcal_data_standard_file = pcal_data_file.split('Pcal')[0] + 'Pcal_standard' + pcal_data_file.split('Pcal')[1]
+ pcal_data_standard = np.loadtxt(pcal_data_standard_file)
+ pcal_data_standard_times = np.transpose(pcal_data_standard)[0]
+ pcal_data = np.loadtxt(pcal_data_file)
+ pcal_data_standard = np.zeros([len(pcal_data_standard), 3])
+ i = 0
+ successes = 0
+ for j in range(len(pcal_data_standard_times)):
+ while pcal_data_standard_times[j] != pcal_data[i][0] and i < len(pcal_data):
+ i += 1
+ if pcal_data_standard_times[j] == pcal_data[i][0]:
+ pcal_data_standard[j] = pcal_data[i]
+ successes += 1
+ if successes != len(pcal_data_standard_times):
+ print("%s: mismatch between data in new file and data in standard file. Expected %d timestamp matches, found %d" % (pcal_data_file, len(pcal_data_standard_times), successes))
+ else:
+ os.system("rm %s" % pcal_data_standard_file)
+ np.savetxt(pcal_data_standard_file, pcal_data_standard)
+
+
+def update_ASD_standard_files():
+ ASD_data_files = os.popen("ls ASD_data/*asd.txt").read().split('\n')[:-1]
+ for ASD_data_file in ASD_data_files:
+ ASD_standard_data_file = ASD_data_file.replace("asd.txt", "asd_standard.txt")
+ ASD_data = np.loadtxt(ASD_data_file)
+ ASD_standard_data = np.loadtxt(ASD_standard_data_file)
+ if len(np.transpose(ASD_data)[0]) != len(np.transpose(ASD_standard_data)[0]):
+ print("%s: mismatch in length of asd standard and new data files (%d != %d)" % (ASD_data_file, len(np.transpose(ASD_standard_data)[0]), len(np.transpose(ASD_data)[0])))
+ elif (np.transpose(ASD_data)[0] != np.transpose(ASD_standard_data)[0]).any()
+ print("%s: frequency vectors don't match" % ASD_data_file)
+ else:
+ os.system("rm %s" % ASD_standard_data_file)
+ np.savetxt(ASD_standard_data_file, ASD_data)
+
+
+def update_TDCF_standard_files():
+ TDCF_data_files = os.popen("ls TDCF_data/*Approx.txt").read().split('\n')[:-1] + os.popen("ls TDCF_data/*Exact.txt").read().split('\n')[:-1]
+ for TDCF_data_file in TDCF_data_files:
+ TDCF_standard_data_file = TDCF_data_files.replace(".txt", "_standard.txt")
+ TDCF_data_standard = np.loadtxt(TDCF_data_standard_file)
+ TDCF_data_standard_times = np.transpose(TDCF_data_standard)[0]
+ TDCF_data = np.loadtxt(TDCF_data_file)
+ TDCF_data_standard = np.zeros([len(TDCF_data_standard), 2])
+ i = 0
+ successes = 0
+ for j in range(len(TDCF_data_standard_times)):
+ while TDCF_data_standard_times[j] != TDCF_data[i][0] and i < len(TDCF_data):
+ i += 1
+ if TDCF_data_standard_times[j] == TDCF_data[i][0]:
+ TDCF_data_standard[j] = TDCF_data[i]
+ successes += 1
+ if successes != len(TDCF_data_standard_times):
+ print("%s: mismatch between data in new file and data in standard file. Expected %d timestamp matches, found %d" % (TDCF_data_file, len(TDCF_data_standard_times), successes))
+ else:
+ os.system("rm %s" % TDCF_data_standard_file)
+ np.savetxt(TDCF_data_standard_file, TDCF_data_standard)
+
+
+def update_STV_standard_files():
+ STV_data_files = ['State_Vector_data/State_Vector_Approx.txt', 'State_Vector_data/State_Vector_Exact.txt']
+ for STV_data_file in STV_data_files:
+ STV_standard_data_file = STV_data_file.replace(".txt", "_standard.txt")
+ STV_data = np.loadtxt(STV_data_file)
+ STV_standard_data = np.loadtxt(STV_standard_data_file)
+ if len(np.transpose(STV_data)[0]) != len(np.transpose(STV_standard_data)[0]):
+ print("%s: mismatch in length of state vector standard and new data files (%d != %d)" % (STV_data_file, len(np.transpose(STV_standard_data)[0]), len(np.transpose(STV_data)[0])))
+ elif (np.transpose(STV_data)[0] != np.transpose(STV_standard_data)[0]).any()
+ print("%s: time vectors don't match" % STV_data_file)
+ else:
+ os.system("rm %s" % STV_standard_data_file)
+ np.savetxt(STV_standard_data_file, STV_data)
+
+