Skip to content
Snippets Groups Projects
Commit ff70b080 authored by Daichi Tsuna's avatar Daichi Tsuna
Browse files

add autocorrelation chi2 calculation 

currently compiles but segfaults in the middle...
parent 536e3fc8
No related branches found
No related tags found
No related merge requests found
Pipeline #50255 passed with warnings
Stage: level0
Stage: level1
Stage: level2
Stage: test-gstlal
Stage: test-inspiral
Stage: test-burst
Stage: test-offline
Hide whitespace changes
Inline Side-by-side
gstlal-burst/bin/gstlal_cs_triggergen
+ 40
− 12
View file @ ff70b080
......@@ -85,6 +85,7 @@ class PSDHandler(simplehandler.Handler):
def __init__(self, mainloop, pipeline, firbank):
simplehandler.Handler.__init__(self,mainloop, pipeline)
self.firbank = firbank
self.triggergen = triggergen
def do_on_message(self, bus, message):
if message.type == Gst.MessageType.ELEMENT and message.get_structure().get_name() == "spectrum":
......@@ -94,57 +95,84 @@ class PSDHandler(simplehandler.Handler):
stability = float(message.src.get_property("n-samples")) / message.src.get_property("average-samples")
if stability > 0.3:
template_bank = [None] * len(template_bank_table)
print >> sys.stderr, "At GPS time", timestamp, "PSD stable"
template_t = [None] * len(template_bank_table)
autocorr = [None] * len(template_bank_table)
# make templates, whiten, put into firbank
for i, row in enumerate(template_bank_table):
# linearly polarized, so just use plus mode time series
template_bank[i], _ = lalsimulation.GenerateStringCusp(1.0,row.central_freq,1.0/options.sample_rate)
template_t[i], _ = lalsimulation.GenerateStringCusp(1.0,row.central_freq,1.0/options.sample_rate)
# zero-pad it to 32 seconds to obtain same deltaF as the PSD
template_bank[i] = lal.ResizeREAL8TimeSeries(template_bank[i],-int(32*options.sample_rate - template_bank[i].data.length)//2,int(32*options.sample_rate))
template_t[i] = lal.ResizeREAL8TimeSeries(template_t[i],-int(32*options.sample_rate - template_t[i].data.length)//2,int(32*options.sample_rate))
# setup of frequency domain
length = template_bank[i].data.length
length = template_t[i].data.length
duration = float(length) / options.sample_rate
epoch = -(length-1) // 2 /options.sample_rate
template_f = lal.CreateCOMPLEX16FrequencySeries("template_freq", LIGOTimeGPS(epoch), psd.f0, 1.0/duration, lal.Unit("strain s"), length // 2 + 1)
fplan = lal.CreateForwardREAL8FFTPlan(length,0)
# FFT to frequency domain
lal.REAL8TimeFreqFFT(template_f,template_bank[i],fplan)
lal.REAL8TimeFreqFFT(template_f,template_t[i],fplan)
# set DC and Nyquist to zero
template_f.data.data[0] = 0.0
template_f.data.data[template_f.data.length-1] = 0.0
# whiten
template_f = lal.WhitenCOMPLEX16FrequencySeries(template_f,psd)
# obtain autocorrelation time series by
# squaring the template and inverse FFTing it
# obtain autocorr time series by squaring template and inverse FFT it
template_f_squared = lal.CreateCOMPLEX16FrequencySeries("whitened template_freq squared", LIGOTimeGPS(epoch), psd.f0, 1.0/duration, lal.Unit("s"), length // 2 + 1)
autocorr[i] = lal.CreateREAL8TimeSeries("autocorr_time", LIGOTimeGPS(epoch), psd.f0, 1.0 / options.sample_rate, lal.Unit("strain"), length)
autocorr_t = lal.CreateREAL8TimeSeries("autocorr_time", LIGOTimeGPS(epoch), psd.f0, 1.0 / options.sample_rate, lal.Unit("strain"), length)
rplan = lal.CreateReverseREAL8FFTPlan(length,0)
template_f_squared.data.data = abs(template_f.data.data)**2
lal.REAL8FreqTimeFFT(autocorr[i],template_f_squared,rplan)
lal.REAL8FreqTimeFFT(autocorr_t,template_f_squared,rplan)
# normalize autocorrelation by central (maximum) value
autocorr[i].data.data /= numpy.max(autocorr[i].data.data)
autocorr_t.data.data /= numpy.max(autocorr_t.data.data)
autocorr_t = autocorr_t.data.data
max_index = numpy.argmax(autocorr_t)
# find the index of the third extremum for the first template
# FIXME we do this because we know that the template with the lowest high-f cutoff will have
# the largest chi2_index.
if i == 0:
extr_ctr = 0
chi2_index = 0
for j in range(max_index+1, len(autocorr_t)):
slope1 = autocorr_t[j+1] - autocorr_t[j]
slope0 = autocorr_t[j] - autocorr_t[j-1]
chi2_index += 1
if(slope1 * slope0 < 0):
extr_ctr += 1
if(extr_ctr == 2):
break
assert extr_ctr == 2, 'could not find 3rd extremum'
# extract the part within the third extremum, setting the peak to be the center.
autocorr[i] = numpy.concatenate((autocorr_t[1:(chi2_index+1)][::-1], autocorr_t[:(chi2_index+1)]))
assert len(autocorr[i])%2==1, 'autocorr must have odd number of samples'
# Inverse FFT template bank back to time domain
template_t[i] = lal.CreateREAL8TimeSeries("whitened template time", LIGOTimeGPS(epoch), psd.f0, 1.0 / options.sample_rate, lal.Unit("strain"), length)
lal.REAL8FreqTimeFFT(template_t[i],template_f,rplan)
# normalize
template_t[i].data.data /= numpy.sqrt(numpy.dot(template_t[i].data.data, template_t[i].data.data))
template_t[i] = template_t[i].data.data
firbank.set_property("latency",-(len(template_t[0]) - 1) // 2)
firbank.set_property("latency", -(len(template_t[0]) - 1) // 2)
firbank.set_property("fir_matrix", template_t)
triggergen.set_property("autocorrelation_matrix", autocorr)
print >> sys.stderr, "finished setting all the properties"
self.firbank = firbank
self.triggergen = triggergen
else:
# use templates with all zeros during burn-in period, that way we won't get any triggers.
print >> sys.stderr, "At GPS time", timestamp, "burn in period"
template = [None] * len(template_bank_table)
autocorr = [None] * len(template_bank_table)
for i, row in enumerate(template_bank_table):
template[i], _ = lalsimulation.GenerateStringCusp(1.0,30,1.0/options.sample_rate)
template[i] = lal.ResizeREAL8TimeSeries(template[i], -int(32*options.sample_rate - template[i].data.length)//2 ,int(32*options.sample_rate))
template[i] = template[i].data.data
template[i] *= 0.0
autocorr[i] = template[i]
firbank.set_property("latency",-(len(template[0]) - 1) // 2)
firbank.set_property("fir_matrix", template)
triggergen.set_property("autocorrelation_matrix", autocorr)
self.firbank = firbank
self.triggergen = triggergen
return True
return False
......@@ -228,7 +256,7 @@ xmldoc.childNodes[-1].appendChild(sngl_burst_table)
# trigger generator
#
head = pipeparts.mkgeneric(pipeline, head, "lal_string_triggergen", threshold = options.threshold, cluster = options.cluster_events, bank_filename = options.template_bank)
head = triggergen = pipeparts.mkgeneric(pipeline, head, "lal_string_triggergen", threshold = options.threshold, cluster = options.cluster_events, bank_filename = options.template_bank)
#
......
gstlal-burst/gst/lal/gstlal_string_triggergen.c
+ 178
− 8
View file @ ff70b080
/*
* Cosmic string trigger generator and autocorrelation chisq plugin element
*/
/*
*======================================================
*
......@@ -24,7 +29,9 @@
#include <glib.h>
#include <gst/gst.h>
#include <gst/audio/audio.h>
#include <gst/base/gstadapter.h>
#include <gst/base/gstbasetransform.h>
#include <gsl/gsl_errno.h>
/*
......@@ -43,9 +50,13 @@
*/
#include <gstlal/gstlal.h>
#include <gstlal_string_triggergen.h>
#include <gstlal/gstaudioadapter.h>
#include <gstlal/gstlal_autocorrelation_chi2.h>
#include <gstlal/gstlal_debug.h>
#include <gstlal/gstlal_peakfinder.h>
/*
*======================================================
......@@ -75,7 +86,7 @@ G_DEFINE_TYPE_WITH_CODE(
*/
#define DEFAULT_THRES 5.5
#define DEFAULT_THRES 4.0
#define DEFAULT_CLUSTER 0.1
......@@ -88,6 +99,19 @@ G_DEFINE_TYPE_WITH_CODE(
*/
static unsigned autocorrelation_length(const GSTLALStringTriggergen *element)
{
return gstlal_autocorrelation_chi2_autocorrelation_length(element->autocorrelation_matrix);
}
static guint64 output_num_bytes(GSTLALStringTriggergen *element)
{
// FIXME don't hardcode sample rate
return (guint64) 8192 * element->adapter->unit_size;
}
static void free_bankfile(GSTLALStringTriggergen *element)
{
g_free(element->bank_filename);
......@@ -146,10 +170,6 @@ static int setup_bankfile_input(GSTLALStringTriggergen *element, char *bank_file
static GstFlowReturn trigger_generator(GSTLALStringTriggergen *element, GstBuffer *inbuf, GstBuffer *outbuf)
{
/*
* find events. we do not use chisq.
*/
GstMapInfo inmap;
float *snrdata;
SnglBurst *triggers = NULL;
......@@ -167,6 +187,17 @@ static GstFlowReturn trigger_generator(GSTLALStringTriggergen *element, GstBuffe
t0 = GST_BUFFER_PTS(inbuf);
length = GST_BUFFER_OFFSET_END(inbuf) - GST_BUFFER_OFFSET(inbuf);
/* copy samples */
gst_audioadapter_copy_samples(element->adapter, element->data, length, NULL, NULL);
/* compute the chisq norm if it doesn't exist */
if (!element->autocorrelation_norm)
element->autocorrelation_norm = gstlal_autocorrelation_chi2_compute_norms(element->autocorrelation_matrix, NULL);
/* check that autocorrelation vector has odd number of samples */
g_assert(autocorrelation_length(element) & 1);
/* find events */
GST_DEBUG_OBJECT(element, "searching %" G_GUINT64_FORMAT " samples at %" GST_TIME_SECONDS_FORMAT " for events with SNR greater than %f", length, GST_TIME_SECONDS_ARGS(t0),element->threshold);
for(sample = 0; sample < length; sample++){
LIGOTimeGPS t;
......@@ -183,6 +214,17 @@ static GstFlowReturn trigger_generator(GSTLALStringTriggergen *element, GstBuffe
*/
element->bank[channel].snr = snr;
element->bank[channel].peak_time = t;
element->bank[channel].chisq_dof = 1.0;
/*
* We calculate chisq each time this update occurs, by defining this as peak.
*/
element->maxdata->values.as_double[channel] = *snrdata;
element->maxdata->samples[channel] = sample;
/* extract data around peak for chisq calculation */
/* put the dat pointer one pad length in */
gstlal_double_series_around_peak(element->maxdata, ((double *) element->data) + element->maxdata->pad * element->num_templates, (double *) element->snr_mat, element->maxdata->pad);
/* calculate chisq */
gstlal_autocorrelation_chi2(&element->bank[channel].chisq, (double complex *) element->snr_mat, autocorrelation_length(element), -((int) autocorrelation_length(element)) / 2, element->threshold, &element->autocorrelation_matrix[channel], NULL, &element->autocorrelation_norm[channel]);
}
} else if(element->bank[channel].snr != 0. && XLALGPSDiff(&t, &element->last_time[channel]) > element->cluster) {
/*
......@@ -192,6 +234,8 @@ static GstFlowReturn trigger_generator(GSTLALStringTriggergen *element, GstBuffe
triggers = g_renew(SnglBurst, triggers, ntriggers + 1);
triggers[ntriggers++] = element->bank[channel];
element->bank[channel].snr = 0.0;
element->bank[channel].chisq = 0.0;
element->bank[channel].chisq_dof = 0.0;
}
}
}
......@@ -291,7 +335,22 @@ static gboolean set_caps(GstBaseTransform *trans, GstCaps *incaps, GstCaps *outc
GSTLALStringTriggergen *element = GSTLAL_STRING_TRIGGERGEN(trans);
gboolean success = gst_audio_info_from_caps(&element->audio_info, incaps);
g_object_set(element->adapter, "unit-size", GST_AUDIO_INFO_WIDTH(&element->audio_info) / 8, NULL);
g_object_set(element->adapter, "unit-size", GST_AUDIO_INFO_WIDTH(&element->audio_info) / 8 * element->num_templates, NULL);
if (element->maxdata)
gstlal_peak_state_free(element->maxdata);
element->maxdata = gstlal_peak_state_new(element->num_templates, GSTLAL_PEAK_DOUBLE_COMPLEX);
/* Update padding any time the autocorrelation property is updated */
if (element->autocorrelation_matrix) {
element->maxdata->pad = autocorrelation_length(element) / 2;
if (element->snr_mat)
free(element->snr_mat);
element->snr_mat = calloc(element->num_templates * autocorrelation_length(element), element->maxdata->unit);
}
/*
* done
*/
return success;
}
......@@ -316,8 +375,16 @@ static gboolean start(GstBaseTransform *trans)
XLALINT8NSToGPS(&element->bank[i].peak_time, 0);
element->bank[i].snr = 0;
/* initialize the last time array, too */
/*
* Initialize the chisq and chisq_dof, too.
* We follow the definition of the previous string search pipeline,
* The actual chi^2 is then chisq/chisq_dof. We can come
* back to the definition later if we have to.
*/
element->bank[i].chisq = 0;
element->bank[i].chisq_dof = 0;
/* initialize the last time array, too */
XLALINT8NSToGPS(&element->last_time[i], 0);
}
......@@ -347,6 +414,17 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
{
GSTLALStringTriggergen *element = GSTLAL_STRING_TRIGGERGEN(trans);
GstFlowReturn result;
guint64 maxsize;
/* The max size to copy from an adapter is the typical output size plus the padding */
maxsize = output_num_bytes(element) + element->adapter->unit_size * element->maxdata->pad * 2;
/* if we haven't allocated storage do it now, we should never try to copy from an adapter with a larger buffer than this */
if (!element->data)
element->data = malloc(maxsize);
/* put the incoming buffer into an adapter */
gst_audioadapter_push(element->adapter, inbuf);
result = trigger_generator(element,inbuf,outbuf);
......@@ -373,7 +451,8 @@ static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuf
enum property {
ARG_THRES = 1,
ARG_CLUSTER,
ARG_BANK_FILENAME
ARG_BANK_FILENAME,
ARG_AUTOCORRELATION_MATRIX
};
......@@ -398,6 +477,33 @@ static void set_property(GObject *object, enum property prop_id, const GValue *v
g_mutex_unlock(&element->bank_lock);
break;
case ARG_AUTOCORRELATION_MATRIX:
g_mutex_lock(&element->bank_lock);
if(element->autocorrelation_matrix)
gsl_matrix_complex_free(element->autocorrelation_matrix);
element->autocorrelation_matrix = gstlal_gsl_matrix_complex_from_g_value_array(g_value_get_boxed(value));
fprintf(stderr,"done setting autocorrelation matrix\n");
/* This should be called any time caps change too */
if(element->maxdata && element->autocorrelation_matrix){
fprintf(stderr, "autocorrelation length = %d\n", autocorrelation_length(element));
element->maxdata->pad = autocorrelation_length(element) / 2;
if (element->snr_mat)
free(element->snr_mat);
element->snr_mat = calloc(element->num_templates * autocorrelation_length(element), element->maxdata->unit);
}
/*
* induce norms to be recomputed
*/
if(element->autocorrelation_norm) {
gsl_vector_free(element->autocorrelation_norm);
element->autocorrelation_norm = NULL;
}
g_mutex_unlock(&element->bank_lock);
break;
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec);
break;
......@@ -427,6 +533,18 @@ static void get_property(GObject *object, enum property prop_id, GValue *value,
g_value_set_string(value, element->bank_filename);
g_mutex_unlock(&element->bank_lock);
break;
case ARG_AUTOCORRELATION_MATRIX:
g_mutex_lock(&element->bank_lock);
if(element->autocorrelation_matrix)
g_value_take_boxed(value, gstlal_g_value_array_from_gsl_matrix_complex(element->autocorrelation_matrix));
else {
GST_WARNING_OBJECT(element, "no autocorrelation matrix");
/* FIXME deprecated.. */
g_value_take_boxed(value, g_value_array_new(0));
}
g_mutex_unlock(&element->bank_lock);
break;
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec);
......@@ -442,12 +560,35 @@ static void finalize(GObject *object)
GSTLALStringTriggergen *element = GSTLAL_STRING_TRIGGERGEN(object);
g_mutex_clear(&element->bank_lock);
free_bankfile(element);
if(element->maxdata) {
gstlal_peak_state_free(element->maxdata);
element->maxdata = NULL;
}
if(element->data){
free(element->data);
element->data = NULL;
}
g_free(element->instrument);
element->instrument = NULL;
g_free(element->channel_name);
element->channel_name = NULL;
gst_audioadapter_clear(element->adapter);
g_object_unref(element->adapter);
if(element->snr_mat) {
free(element->snr_mat);
element->snr_mat = NULL;
}
if(element->autocorrelation_matrix) {
gsl_matrix_complex_free(element->autocorrelation_matrix);
element->autocorrelation_matrix = NULL;
}
if(element->autocorrelation_norm) {
gsl_vector_free(element->autocorrelation_norm);
element->autocorrelation_norm = NULL;
}
G_OBJECT_CLASS(gstlal_string_triggergen_parent_class)->finalize(object);
}
......@@ -538,6 +679,30 @@ static void gstlal_string_triggergen_class_init(GSTLALStringTriggergenClass *kla
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
)
);
g_object_class_install_property(
gobject_class,
ARG_AUTOCORRELATION_MATRIX,
g_param_spec_value_array(
"autocorrelation-matrix",
"Autocorrelation Matrix",
"Array of autocorrelation vectors. Number of vectors (rows) in matrix sets number of channels. All vectors must have the same length.",
g_param_spec_value_array(
"autocorrelation",
"Autocorrelation",
"Array of autocorrelation samples.",
/* FIXME: should be complex */
g_param_spec_double(
"sample",
"Sample",
"Autocorrelation sample",
-G_MAXDOUBLE, G_MAXDOUBLE, 0.0,
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS
),
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS
),
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS
)
);
}
/*
......@@ -551,10 +716,15 @@ static void gstlal_string_triggergen_init(GSTLALStringTriggergen *element)
element->adapter = g_object_new(GST_TYPE_AUDIOADAPTER, NULL);
element->bank_filename = NULL;
element->bank = NULL;
element->data = NULL;
element->maxdata = NULL;
element->instrument = NULL;
element->channel_name = NULL;
element->num_templates = 0;
element->last_time = NULL;
element->snr_mat = NULL;
element->audio_info.bpf = 0; /* impossible value */
element->autocorrelation_matrix = NULL;
element->autocorrelation_norm = NULL;
gst_base_transform_set_gap_aware(GST_BASE_TRANSFORM(element), TRUE);
}
gstlal-burst/gst/lal/gstlal_string_triggergen.h
+ 11
− 5
View file @ ff70b080
......@@ -5,9 +5,14 @@
#include <glib.h>
#include <gst/gst.h>
#include <gst/audio/audio.h>
#include <gst/base/gstadapter.h>
#include <gst/base/gstbasetransform.h>
#include <gstlal/gstaudioadapter.h>
#include <gstlal/gstlal_peakfinder.h>
#include <lal/LIGOMetadataTables.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_matrix_float.h>
G_BEGIN_DECLS
......@@ -40,20 +45,21 @@ typedef struct {
GstAudioInfo audio_info;
/*
* extracting triggers above threshold
*/
float threshold;
float cluster;
GMutex bank_lock;
gsl_matrix_complex *autocorrelation_matrix;
gsl_vector *autocorrelation_norm;
char *bank_filename;
SnglBurst *bank;
void *data;
struct gstlal_peak_state *maxdata;
gchar *instrument;
gchar *channel_name;
SnglBurst *bank;
gint num_templates;
LIGOTimeGPS *last_time;
void *snr_mat;
} GSTLALStringTriggergen;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment