From e20a62465e0594ba0bcf548525a88ae600d3fe57 Mon Sep 17 00:00:00 2001
From: Patrick Godwin <patrick.godwin@ligo.org>
Date: Sun, 26 Aug 2018 09:39:24 -0700
Subject: [PATCH] gstlal_sngltrigger: factor out code that adds single peaks to
buffers to avoid looping through channels to select max snr channel when
using this option
---
gstlal-burst/lib/gstlal_sngltrigger.c | 218 ++++++++++++++------------
gstlal-burst/lib/gstlal_sngltrigger.h | 1 +
2 files changed, 115 insertions(+), 104 deletions(-)
diff --git a/gstlal-burst/lib/gstlal_sngltrigger.c b/gstlal-burst/lib/gstlal_sngltrigger.c
index 9cb1dc32f6..30547d3047 100644
--- a/gstlal-burst/lib/gstlal_sngltrigger.c
+++ b/gstlal-burst/lib/gstlal_sngltrigger.c
@@ -1,5 +1,6 @@
/*
* Copyright (C) 2012,2013,2015,2016 Chad Hanna
+ * 2017-2018 Duncan Meacher
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
@@ -57,120 +58,129 @@ GstBuffer *gstlal_sngltrigger_new_buffer_from_peak(struct gstlal_peak_state *inp
GST_BUFFER_DURATION(srcbuf) = (GstClockTime) gst_util_uint64_scale_int_round(GST_SECOND, length, rate);
if (input->num_events) {
- if(max_snr)
- {
- max_snr_index = gstlal_peak_max_over_channels(input);
- /* Type casting unsigned int (guint) to int */
- g_assert(max_snr_index >= 0 && max_snr_index < (int)input->channels);
- }
- guint channel;
- for(channel = 0; channel < input->channels; channel++) {
- struct GSTLALSnglTrigger *event;
- SnglTriggerTable *parent;
- double complex maxdata_channel = 0;
+ /* create a single buffer from max snr channel */
+ if(max_snr) {
+ max_snr_index = gstlal_peak_max_over_channels(input);
/* Type casting unsigned int (guint) to int */
- if(max_snr && (int)channel != max_snr_index)
- continue;
+ g_assert(max_snr_index >= 0 && max_snr_index < (int)input->channels);
- switch (input->type)
- {
- case GSTLAL_PEAK_COMPLEX:
- maxdata_channel = (double complex) input->interpvalues.as_float_complex[channel];
- break;
+ /* create a single buffer from max snr channel */
+ add_buffer_from_channel(input, channel_name, pad, length, time, rate, chi2, snr_matrix_view, max_snr_index, srcbuf);
- case GSTLAL_PEAK_DOUBLE_COMPLEX:
- maxdata_channel = (double complex) input->interpvalues.as_double_complex[channel];
- break;
+ /* add a buffer from each channel */
+ } else {
- default:
- g_assert(input->type == GSTLAL_PEAK_COMPLEX || input->type == GSTLAL_PEAK_DOUBLE_COMPLEX);
+ guint channel;
+ for(channel = 0; channel < input->channels; channel++) {
+ add_buffer_from_channel(input, channel_name, pad, length, time, rate, chi2, snr_matrix_view, channel, srcbuf);
}
+ }
+ }
- if (!maxdata_channel)
- continue;
-
- /*
- * allocate new event structure
- */
-
-
- /*
- * Populate the SNR snippet if available
- * FIXME: only supported for single precision at the moment
- */
- if (snr_matrix_view)
- {
- /* Get the column of SNR we are interested in */
- gsl_vector_complex_float_view snr_vector_view = gsl_matrix_complex_float_column(&(snr_matrix_view->matrix), channel);
- /* Allocate an empty time series to hold it. The event takes ownership, so no need to free it*/
- event = gstlal_sngltrigger_new(snr_vector_view.vector.size);
- /* Make a GSL view of the time series array data */
- gsl_vector_complex_float_view snr_series_view = gsl_vector_complex_float_view_array((float *) event->snr, event->length);
- /* Use BLAS to do the copy */
- gsl_blas_ccopy (&(snr_vector_view.vector), &(snr_series_view.vector));
- } else
- event = gstlal_sngltrigger_new(0);
-
- parent = (SnglTriggerTable *) event;
- if (!event) {
- /* FIXME handle error */
- }
+ return srcbuf;
+}
- /*
- * populate
- */
-
- //*parent = bankarray[channel];
- strcpy(parent->channel, channel_name);
- parent->snr = cabs(maxdata_channel);
- parent->phase = carg(maxdata_channel);
- parent->channel_index = channel;
-
- XLALINT8NSToGPS(&event->epoch, time);
- {
- LIGOTimeGPS end_time = event->epoch;
- XLALGPSAdd(&end_time, (double) input->interpsamples[channel] / rate);
- XLALGPSAddGPS(&parent->end, &end_time);
- }
- XLALGPSAdd(&event->epoch, (double) (input->samples[channel] - input->pad) / rate);
- event->deltaT = 1. / rate;
-
- /* populate chi squared if we have it */
- parent->chisq = 0.0;
- switch (input->type)
- {
- case GSTLAL_PEAK_COMPLEX:
- if (chi2) parent->chisq = (double) *(((float *) chi2 ) + channel);
- break;
-
- case GSTLAL_PEAK_DOUBLE_COMPLEX:
- if (chi2) parent->chisq = (double) *(((double *) chi2 ) + channel);
- break;
-
- default:
- g_assert(input->type == GSTLAL_PEAK_COMPLEX || input->type == GSTLAL_PEAK_DOUBLE_COMPLEX);
- }
+void add_buffer_from_channel(struct gstlal_peak_state *input, char *channel_name, GstPad *pad, guint64 length, GstClockTime time, guint rate, void *chi2, gsl_matrix_complex_float_view *snr_matrix_view, int channel, GstBuffer *srcbuf)
+{
+ struct GSTLALSnglTrigger *event;
+ SnglTriggerTable *parent;
+ double complex maxdata_channel = 0;
+
+ switch (input->type)
+ {
+ case GSTLAL_PEAK_COMPLEX:
+ maxdata_channel = (double complex) input->interpvalues.as_float_complex[channel];
+ break;
+
+ case GSTLAL_PEAK_DOUBLE_COMPLEX:
+ maxdata_channel = (double complex) input->interpvalues.as_double_complex[channel];
+ break;
+
+ default:
+ g_assert(input->type == GSTLAL_PEAK_COMPLEX || input->type == GSTLAL_PEAK_DOUBLE_COMPLEX);
+ }
- /*
- * add to buffer
- */
-
- gst_buffer_append_memory(
- srcbuf,
- gst_memory_new_wrapped(
- GST_MEMORY_FLAG_READONLY | GST_MEMORY_FLAG_PHYSICALLY_CONTIGUOUS,
- event,
- sizeof(*event) + sizeof(event->snr[0]) * event->length,
- 0,
- sizeof(*event) + sizeof(event->snr[0]) * event->length,
- event,
- (GDestroyNotify) gstlal_sngltrigger_free
- )
- );
- }
+ if (!maxdata_channel)
+ return;
+
+ /*
+ * allocate new event structure
+ */
+
+
+ /*
+ * Populate the SNR snippet if available
+ * FIXME: only supported for single precision at the moment
+ */
+ if (snr_matrix_view)
+ {
+ /* Get the column of SNR we are interested in */
+ gsl_vector_complex_float_view snr_vector_view = gsl_matrix_complex_float_column(&(snr_matrix_view->matrix), channel);
+ /* Allocate an empty time series to hold it. The event takes ownership, so no need to free it*/
+ event = gstlal_sngltrigger_new(snr_vector_view.vector.size);
+ /* Make a GSL view of the time series array data */
+ gsl_vector_complex_float_view snr_series_view = gsl_vector_complex_float_view_array((float *) event->snr, event->length);
+ /* Use BLAS to do the copy */
+ gsl_blas_ccopy (&(snr_vector_view.vector), &(snr_series_view.vector));
+ } else
+ event = gstlal_sngltrigger_new(0);
+
+ parent = (SnglTriggerTable *) event;
+ if (!event) {
+ /* FIXME handle error */
}
- return srcbuf;
+ /*
+ * populate
+ */
+
+ //*parent = bankarray[channel];
+ strcpy(parent->channel, channel_name);
+ parent->snr = cabs(maxdata_channel);
+ parent->phase = carg(maxdata_channel);
+ parent->channel_index = channel;
+
+ XLALINT8NSToGPS(&event->epoch, time);
+ {
+ LIGOTimeGPS end_time = event->epoch;
+ XLALGPSAdd(&end_time, (double) input->interpsamples[channel] / rate);
+ XLALGPSAddGPS(&parent->end, &end_time);
+ }
+ XLALGPSAdd(&event->epoch, (double) (input->samples[channel] - input->pad) / rate);
+ event->deltaT = 1. / rate;
+
+ /* populate chi squared if we have it */
+ parent->chisq = 0.0;
+ switch (input->type)
+ {
+ case GSTLAL_PEAK_COMPLEX:
+ if (chi2) parent->chisq = (double) *(((float *) chi2 ) + channel);
+ break;
+
+ case GSTLAL_PEAK_DOUBLE_COMPLEX:
+ if (chi2) parent->chisq = (double) *(((double *) chi2 ) + channel);
+ break;
+
+ default:
+ g_assert(input->type == GSTLAL_PEAK_COMPLEX || input->type == GSTLAL_PEAK_DOUBLE_COMPLEX);
+ }
+
+ /*
+ * add to buffer
+ */
+
+ gst_buffer_append_memory(
+ srcbuf,
+ gst_memory_new_wrapped(
+ GST_MEMORY_FLAG_READONLY | GST_MEMORY_FLAG_PHYSICALLY_CONTIGUOUS,
+ event,
+ sizeof(*event) + sizeof(event->snr[0]) * event->length,
+ 0,
+ sizeof(*event) + sizeof(event->snr[0]) * event->length,
+ event,
+ (GDestroyNotify) gstlal_sngltrigger_free
+ )
+ );
}
diff --git a/gstlal-burst/lib/gstlal_sngltrigger.h b/gstlal-burst/lib/gstlal_sngltrigger.h
index f3ba91e256..04a4cac5b1 100644
--- a/gstlal-burst/lib/gstlal_sngltrigger.h
+++ b/gstlal-burst/lib/gstlal_sngltrigger.h
@@ -65,6 +65,7 @@ SnglTriggerTable;
*/
GstBuffer *gstlal_sngltrigger_new_buffer_from_peak(struct gstlal_peak_state *input, char *channel_name, GstPad *pad, guint64 offset, guint64 length, GstClockTime time, guint rate, void *chi2, gsl_matrix_complex_float_view *snr_matrix_view, GstClockTimeDiff, gboolean max_snr);
+void add_buffer_from_channel(struct gstlal_peak_state *input, char *channel_name, GstPad *pad, guint64 length, GstClockTime time, guint rate, void *chi2, gsl_matrix_complex_float_view *snr_matrix_view, int channel, GstBuffer *srcbuf);
G_END_DECLS
#endif /* __GSTLAL_SNGLTRIGGER_H__ */
--
GitLab