Commit 1ae71845 authored by Ryan Magee's avatar Ryan Magee
Browse files

Bankchisq: This commit incorporates the bankchisq calculation in itacac and...

Bankchisq: This commit incorporates the bankchisq calculation in itacac and modifies python modules that used the column to store intermediate values. This commit does NOT include necessary changes to store the bankchisq histograms, and it does not incorporate it into the ranking statistic
parent b1ee4e8c
Pipeline #141297 passed with stages
in 28 minutes and 50 seconds
......@@ -278,6 +278,9 @@ back += report.Description("The chi-squared test checks that the snr accumulated
imgfooter = "Chi-squared vs snr for single detectors after coincidence. Blue points are full data zero lag, red are software injections and black are time slides."
back += report.ImageGrid("Chi-squared Distributions", grid_size=3, footer=imgfooter).glob(os.path.join(opts.glob_path, '*3_chi2_vs_rho_*closedbox*.png'))
imgfooter = "Bank chi-squared vs snr for single detectors after coincidence. Blue points are full data zero lag, red are software injections and black are time slides."
back += report.ImageGrid("Bank chi-squared Distributions", grid_size=3, footer=imgfooter).glob(os.path.join(opts.glob_path, '*3_bankchi2_vs_rho_*closedbox*.png'))
imgfooter = "Comparison of SNR in pairs of detectors."
back += report.ImageGrid("Signal-to-Noise Ratio", grid_size=3, footer=imgfooter).glob(os.path.join(opts.glob_path, '*4_rho_*_vs_*closedbox*.png'))
......@@ -290,6 +293,9 @@ if opts.open_box:
imgfooter = "Chi-squared vs snr for single detectors after coincidence. Blue points are full data zero lag, red are software injections and black are time slides."
back += report.ImageGrid("Chi-squared Distributions: Zero lag", grid_size=3, footer=imgfooter).glob(os.path.join(opts.glob_path, '*3_chi2_vs_rho_*openbox*.png'))
imgfooter = "Bank chi-squared vs snr for single detectors after coincidence. Blue points are full data zero lag, red are software injections and black are time slides."
back += report.ImageGrid("Bank chi-squared Distributions: Zero lag", grid_size=3, footer=imgfooter).glob(os.path.join(opts.glob_path, '*3_bankchi2_vs_rho_*openbox*.png'))
imgfooter = "Comparison of SNR in pairs of detectors."
back += report.ImageGrid("Signal-to-Noise Ratio: Zero lag", grid_size=3, footer=imgfooter).glob(os.path.join(opts.glob_path, '*4_rho_*_vs_*openbox*.png'))
......@@ -396,8 +402,12 @@ summary += report.Description("The chi-squared test checks that the snr accumula
imgfooter = "Chi-squared vs snr for single detectors after coincidence. Blue points are full data zero lag, red are software injections and black are time slides."
summary += report.ImageGrid("Chi-squared Distributions", grid_size=3, footer=imgfooter).glob(os.path.join(opts.glob_path, '*3_chi2_vs_rho_*closedbox*.png'))
imgfooter = "Bank chi-squared vs snr for single detectors after coincidence. Blue points are full data zero lag, red are software injections and black are time slides."
summary += report.ImageGrid("Bank chi-squared Distributions", grid_size=3, footer=imgfooter).glob(os.path.join(opts.glob_path, '*3_bankchi2_vs_rho_*closedbox*.png'))
if opts.open_box:
summary += report.ImageGrid("Chi-squared Distributions: Zero lag", grid_size=3, footer=imgfooter).glob(os.path.join(opts.glob_path, '*3_chi2_vs_rho_*openbox*.png'))
summary += report.ImageGrid("Bank chi-squared Distributions: Zero lag", grid_size=3, footer=imgfooter).glob(os.path.join(opts.glob_path, '*3_bankchi2_vs_rho_*openbox*.png'))
## Include IFAR, lnL plots + summary table
summary += report.Header('Money Plots')
......
This diff is collapsed.
......@@ -88,6 +88,8 @@ typedef struct {
struct data_container *data;
void *chi2;
void *tmp_chi2;
void *bankchi2;
void *tmp_bankchi2;
char *bank_filename;
char *instrument;
char *channel_name;
......@@ -99,6 +101,8 @@ typedef struct {
gsl_matrix_complex *autocorrelation_matrix;
gsl_matrix_int *autocorrelation_mask;
gsl_vector *autocorrelation_norm;
gsl_matrix_complex *bankcorrelation_matrix;
gsl_vector *bankcorrelation_norm;
void *snr_mat;
void *tmp_snr_mat;
gsl_matrix_complex_float_view snr_matrix_view;
......
......@@ -328,7 +328,7 @@ int gstlal_set_min_offset_in_snglinspiral_array(SnglInspiralTable *bankarray, in
return 0;
}
int populate_snglinspiral_buffer(GstBuffer *srcbuf, struct gstlal_peak_state *input, SnglInspiralTable *bankarray, GstPad *pad, guint64 length, GstClockTime time, guint rate, void *chi2, gsl_matrix_complex_float_view *G1_snr_matrix_view, gsl_matrix_complex_float_view *H1_snr_matrix_view, gsl_matrix_complex_float_view *K1_snr_matrix_view, gsl_matrix_complex_float_view *L1_snr_matrix_view, gsl_matrix_complex_float_view *V1_snr_matrix_view)
int populate_snglinspiral_buffer(GstBuffer *srcbuf, struct gstlal_peak_state *input, SnglInspiralTable *bankarray, GstPad *pad, guint64 length, GstClockTime time, guint rate, void *chi2, void *bankchi2, gsl_matrix_complex_float_view *G1_snr_matrix_view, gsl_matrix_complex_float_view *H1_snr_matrix_view, gsl_matrix_complex_float_view *K1_snr_matrix_view, gsl_matrix_complex_float_view *L1_snr_matrix_view, gsl_matrix_complex_float_view *V1_snr_matrix_view)
{
guint channel;
guint G1_snr_timeseries_length, H1_snr_timeseries_length, K1_snr_timeseries_length, L1_snr_timeseries_length, V1_snr_timeseries_length;
......@@ -455,14 +455,18 @@ int populate_snglinspiral_buffer(GstBuffer *srcbuf, struct gstlal_peak_state *in
/* populate chi squared if we have it */
parent->chisq = 0.0;
parent->chisq_dof = 1;
parent->bank_chisq = 0.0;
parent->bank_chisq_dof = 1;
switch (input->type)
{
case GSTLAL_PEAK_COMPLEX:
if (chi2) parent->chisq = (double) *(((float *) chi2 ) + channel);
if (bankchi2) parent->bank_chisq = (double) *(((float *) bankchi2 ) + channel);
break;
case GSTLAL_PEAK_DOUBLE_COMPLEX:
if (chi2) parent->chisq = (double) *(((double *) chi2 ) + channel);
if (bankchi2) parent->bank_chisq = (double) *(((double *) bankchi2 ) + channel);
break;
default:
......@@ -490,7 +494,7 @@ int populate_snglinspiral_buffer(GstBuffer *srcbuf, struct gstlal_peak_state *in
return 0;
}
GstBuffer *gstlal_snglinspiral_new_buffer_from_peak(struct gstlal_peak_state *input, SnglInspiralTable *bankarray, GstPad *pad, guint64 offset, guint64 length, GstClockTime time, guint rate, void *chi2, gsl_matrix_complex_float_view *G1_snr_matrix_view, gsl_matrix_complex_float_view *H1_snr_matrix_view, gsl_matrix_complex_float_view *K1_snr_matrix_view, gsl_matrix_complex_float_view *L1_snr_matrix_view, gsl_matrix_complex_float_view *V1_snr_matrix_view, GstClockTimeDiff timediff)
GstBuffer *gstlal_snglinspiral_new_buffer_from_peak(struct gstlal_peak_state *input, SnglInspiralTable *bankarray, GstPad *pad, guint64 offset, guint64 length, GstClockTime time, guint rate, void *chi2, void *bankchi2, gsl_matrix_complex_float_view *G1_snr_matrix_view, gsl_matrix_complex_float_view *H1_snr_matrix_view, gsl_matrix_complex_float_view *K1_snr_matrix_view, gsl_matrix_complex_float_view *L1_snr_matrix_view, gsl_matrix_complex_float_view *V1_snr_matrix_view, GstClockTimeDiff timediff)
{
GstBuffer *srcbuf = gst_buffer_new();
......@@ -512,12 +516,12 @@ GstBuffer *gstlal_snglinspiral_new_buffer_from_peak(struct gstlal_peak_state *in
GST_BUFFER_DURATION(srcbuf) = (GstClockTime) gst_util_uint64_scale_int_round(GST_SECOND, length, rate);
if (input->num_events || input->no_peaks_past_threshold) {
populate_snglinspiral_buffer(srcbuf, input, bankarray, pad, length, time, rate, chi2, G1_snr_matrix_view, H1_snr_matrix_view, K1_snr_matrix_view, L1_snr_matrix_view, V1_snr_matrix_view);
populate_snglinspiral_buffer(srcbuf, input, bankarray, pad, length, time, rate, chi2, bankchi2, G1_snr_matrix_view, H1_snr_matrix_view, K1_snr_matrix_view, L1_snr_matrix_view, V1_snr_matrix_view);
}
return srcbuf;
}
int gstlal_snglinspiral_append_peak_to_buffer(GstBuffer *srcbuf, struct gstlal_peak_state *input, SnglInspiralTable *bankarray, GstPad *pad, guint64 offset, guint64 length, GstClockTime time, guint rate, void *chi2, gsl_matrix_complex_float_view *G1_snr_matrix_view, gsl_matrix_complex_float_view *H1_snr_matrix_view, gsl_matrix_complex_float_view *K1_snr_matrix_view, gsl_matrix_complex_float_view *L1_snr_matrix_view, gsl_matrix_complex_float_view *V1_snr_matrix_view)
int gstlal_snglinspiral_append_peak_to_buffer(GstBuffer *srcbuf, struct gstlal_peak_state *input, SnglInspiralTable *bankarray, GstPad *pad, guint64 offset, guint64 length, GstClockTime time, guint rate, void *chi2, void *bankchi2, gsl_matrix_complex_float_view *G1_snr_matrix_view, gsl_matrix_complex_float_view *H1_snr_matrix_view, gsl_matrix_complex_float_view *K1_snr_matrix_view, gsl_matrix_complex_float_view *L1_snr_matrix_view, gsl_matrix_complex_float_view *V1_snr_matrix_view)
{
//
// Add peak information to a buffer, GST_BUFFER_OFFSET cannot be
......@@ -530,7 +534,7 @@ int gstlal_snglinspiral_append_peak_to_buffer(GstBuffer *srcbuf, struct gstlal_p
GST_BUFFER_DURATION(srcbuf) = (GstClockTime) gst_util_uint64_scale_int_round(GST_SECOND, GST_BUFFER_OFFSET_END(srcbuf) - GST_BUFFER_OFFSET(srcbuf), rate);
}
populate_snglinspiral_buffer(srcbuf, input, bankarray, pad, length, time, rate, chi2, G1_snr_matrix_view, H1_snr_matrix_view, K1_snr_matrix_view, L1_snr_matrix_view, V1_snr_matrix_view);
populate_snglinspiral_buffer(srcbuf, input, bankarray, pad, length, time, rate, chi2, bankchi2, G1_snr_matrix_view, H1_snr_matrix_view, K1_snr_matrix_view, L1_snr_matrix_view, V1_snr_matrix_view);
return 0;
}
......@@ -55,8 +55,8 @@ void gstlal_snglinspiral_array_free(SnglInspiralTable *bankarray);
/*
* FIXME: only support single precision SNR snippets at the moment
*/
GstBuffer *gstlal_snglinspiral_new_buffer_from_peak(struct gstlal_peak_state *input, SnglInspiralTable *bankarray, GstPad *pad, guint64 offset, guint64 length, GstClockTime time, guint rate, void *chi2, gsl_matrix_complex_float_view *G1_snr_matrix_view, gsl_matrix_complex_float_view *H1_snr_matrix_view, gsl_matrix_complex_float_view *K1_snr_matrix_view, gsl_matrix_complex_float_view *L1_snr_matrix_view, gsl_matrix_complex_float_view *V1_snr_matrix_view, GstClockTimeDiff);
int gstlal_snglinspiral_append_peak_to_buffer(GstBuffer *srcbuf, struct gstlal_peak_state *input, SnglInspiralTable *bankarray, GstPad *pad, guint64 offset, guint64 length, GstClockTime time, guint rate, void *chi2, gsl_matrix_complex_float_view *G1_snr_matrix_view, gsl_matrix_complex_float_view *H1_snr_matrix_view, gsl_matrix_complex_float_view *K1_snr_matrix_view, gsl_matrix_complex_float_view *L1_snr_matrix_view, gsl_matrix_complex_float_view *V1_snr_matrix_view);
GstBuffer *gstlal_snglinspiral_new_buffer_from_peak(struct gstlal_peak_state *input, SnglInspiralTable *bankarray, GstPad *pad, guint64 offset, guint64 length, GstClockTime time, guint rate, void *chi2, void *bankchi2, gsl_matrix_complex_float_view *G1_snr_matrix_view, gsl_matrix_complex_float_view *H1_snr_matrix_view, gsl_matrix_complex_float_view *K1_snr_matrix_view, gsl_matrix_complex_float_view *L1_snr_matrix_view, gsl_matrix_complex_float_view *V1_snr_matrix_view, GstClockTimeDiff);
int gstlal_snglinspiral_append_peak_to_buffer(GstBuffer *srcbuf, struct gstlal_peak_state *input, SnglInspiralTable *bankarray, GstPad *pad, guint64 offset, guint64 length, GstClockTime time, guint rate, void *chi2, void *bankchi2, gsl_matrix_complex_float_view *G1_snr_matrix_view, gsl_matrix_complex_float_view *H1_snr_matrix_view, gsl_matrix_complex_float_view *K1_snr_matrix_view, gsl_matrix_complex_float_view *L1_snr_matrix_view, gsl_matrix_complex_float_view *V1_snr_matrix_view);
G_END_DECLS
......
......@@ -312,7 +312,7 @@ def chisq_distribution(df, non_centralities, size):
class CoincsDocument(object):
sngl_inspiral_columns = ("process:process_id", "ifo", "end_time", "end_time_ns", "eff_distance", "coa_phase", "mass1", "mass2", "snr", "chisq", "chisq_dof", "bank_chisq", "bank_chisq_dof", "sigmasq", "spin1x", "spin1y", "spin1z", "spin2x", "spin2y", "spin2z", "template_duration", "event_id", "Gamma0", "Gamma1")
sngl_inspiral_columns = ("process:process_id", "ifo", "end_time", "end_time_ns", "eff_distance", "coa_phase", "mass1", "mass2", "snr", "chisq", "chisq_dof", "bank_chisq", "bank_chisq_dof", "sigmasq", "spin1x", "spin1y", "spin1z", "spin2x", "spin2y", "spin2z", "template_duration", "event_id", "Gamma0", "Gamma1", "Gamma2")
def __init__(self, url, process_params, process_start_time, comment, instruments, seg, offsetvectors, injection_filename = None, tmp_path = None, replace_file = None, verbose = False):
#
......
......@@ -993,9 +993,9 @@ class Handler(simplehandler.Handler):
# FIXME: ugly way to get the instrument
instruments = set([event.ifo for event in events])
# FIXME calculate a chisq weighted SNR and store it in the bank_chisq column
# FIXME calculate a chisq weighted SNR and store it in the Gamma9 column
for event in events:
event.bank_chisq = event.snr / ((1 + max(1., event.chisq)**3)/2.0)**(1./5.)
event.Gamma2 = event.snr / ((1 + max(1., event.chisq)**3)/2.0)**(1./5.)
# extract segment. move the segment's upper
# boundary to include all triggers. ARGH the 1 ns
......
......@@ -702,10 +702,10 @@ def mkLLOIDmulti(pipeline, detectors, banks, psd, psd_fft_length = 32, ht_gate_t
head = itacac_dict[bank.bank_id]
pad = head.get_request_pad("sink%d" % len(head.sinkpads))
if instrument == 'H1' or instrument == 'L1':
for prop, val in [("n", nsamps_window), ("snr-thresh", LnLRDensity.snr_min), ("bank_filename", bank.template_bank_filename), ("sigmasq", bank.sigmasq), ("autocorrelation_matrix", pipeio.repack_complex_array_to_real(bank.autocorrelation_bank)), ("autocorrelation_mask", bank.autocorrelation_mask)]:
for prop, val in [("n", nsamps_window), ("snr-thresh", LnLRDensity.snr_min), ("bank_filename", bank.template_bank_filename), ("sigmasq", bank.sigmasq), ("autocorrelation_matrix", pipeio.repack_complex_array_to_real(bank.autocorrelation_bank)), ("autocorrelation_mask", bank.autocorrelation_mask), ("bankcorrelation_matrix", pipeio.repack_complex_array_to_real(bank.bank_correlation_matrix))]:
pad.set_property(prop, val)
else:
for prop, val in [("n", nsamps_window), ("snr-thresh", LnLRDensity.snr_min), ("bank_filename", bank.template_bank_filename), ("sigmasq", bank.sigmasq), ("autocorrelation_matrix", pipeio.repack_complex_array_to_real(bank.autocorrelation_bank)), ("autocorrelation_mask", bank.autocorrelation_mask)]:
for prop, val in [("n", nsamps_window), ("snr-thresh", LnLRDensity.snr_min), ("bank_filename", bank.template_bank_filename), ("sigmasq", bank.sigmasq), ("autocorrelation_matrix", pipeio.repack_complex_array_to_real(bank.autocorrelation_bank)), ("autocorrelation_mask", bank.autocorrelation_mask), ("bankcorrelation_matrix", pipeio.repack_complex_array_to_real(bank.bank_correlation_matrix))]:
pad.set_property(prop, val)
else:
raise NotImplementedError("Currently only 'autochisq' is supported")
......
......@@ -184,10 +184,19 @@ class Bank(object):
# Assign template banks to fragments
self.bank_fragments = [BankFragment(rate,begin,end) for rate,begin,end in time_slices]
# Setup bank correlation matrix
self.bank_correlation_matrix = None
for i, bank_fragment in enumerate(self.bank_fragments):
if verbose:
print >>sys.stderr, "constructing template decomposition %d of %d: %g s ... %g s" % (i + 1, len(self.bank_fragments), -bank_fragment.end, -bank_fragment.start)
bank_fragment.set_template_bank(template_bank[i], tolerance, self.snr_threshold, identity_transform = identity_transform, verbose = verbose)
cmix = bank_fragment.mix_matrix[:,::2] + 1.j * bank_fragment.mix_matrix[:,1::2]
if self.bank_correlation_matrix is None:
self.bank_correlation_matrix = numpy.dot(numpy.conj(cmix.T), cmix)
else:
self.bank_correlation_matrix += numpy.dot(numpy.conj(cmix.T), cmix)
if bank_fragment.sum_of_squares_weights is not None:
self.gate_threshold = sum_of_squares_threshold_from_fap(gate_fap, numpy.array([weight**2 for bank_fragment in self.bank_fragments for weight in bank_fragment.sum_of_squares_weights], dtype = "double"))
......@@ -317,6 +326,7 @@ def write_bank(filename, banks, psd_input, cliplefts = None, cliprights = None,
bank.autocorrelation_bank = bank.autocorrelation_bank[clipleft:clipright,:]
bank.autocorrelation_mask = bank.autocorrelation_mask[clipleft:clipright,:]
bank.sigmasq = bank.sigmasq[clipleft:clipright]
bank.bank_correlation_matrix = bank.bank_correlation_matrix[clipleft:clipright,clipleft:clipright]
# Add root-level arrays
# FIXME: ligolw format now supports complex-valued data
......@@ -324,6 +334,8 @@ def write_bank(filename, banks, psd_input, cliplefts = None, cliprights = None,
root.appendChild(ligolw_array.Array.build('autocorrelation_bank_imag', bank.autocorrelation_bank.imag))
root.appendChild(ligolw_array.Array.build('autocorrelation_mask', bank.autocorrelation_mask))
root.appendChild(ligolw_array.Array.build('sigmasq', numpy.array(bank.sigmasq)))
root.appendChild(ligolw_array.Array.build('bank_correlation_matrix_real', bank.bank_correlation_matrix.real))
root.appendChild(ligolw_array.Array.build('bank_correlation_matrix_imag', bank.bank_correlation_matrix.imag))
# Write bank fragments
for i, frag in enumerate(bank.bank_fragments):
......@@ -407,6 +419,9 @@ def read_banks(filename, contenthandler, verbose = False):
bank.autocorrelation_bank = ligolw_array.get_array(root, 'autocorrelation_bank_real').array + 1j * ligolw_array.get_array(root, 'autocorrelation_bank_imag').array
bank.autocorrelation_mask = ligolw_array.get_array(root, 'autocorrelation_mask').array
bank.sigmasq = ligolw_array.get_array(root, 'sigmasq').array
bank_correlation_real = ligolw_array.get_array(root, 'bank_correlation_matrix_real').array
bank_correlation_imag = ligolw_array.get_array(root, 'bank_correlation_matrix_imag').array
bank.bank_correlation_matrix = bank_correlation_real + 1j * bank_correlation_imag
# prepare the horizon distance factors
bank.horizon_factors = dict((row.template_id, sigmasq**.5) for row, sigmasq in zip(bank.sngl_inspiral_table, bank.sigmasq))
......
......@@ -184,7 +184,7 @@ CREATE TEMPORARY TABLE _cluster_info_ AS
(coinc_inspiral.end_time - (SELECT MIN(end_time) FROM coinc_inspiral)) + 1e-9 * coinc_inspiral.end_time_ns AS end_time,
coinc_inspiral.snr AS snr,
-- NOTE FIXME we are using bank_chisq to hold a chisq weighted snr value
(SELECT SUM(sngl.bank_chisq * sngl.bank_chisq) FROM sngl_inspiral as sngl WHERE sngl.event_id IN (SELECT event_id FROM coinc_event_map WHERE coinc_event_map.coinc_event_id == coinc_event.coinc_event_id)) AS ranking_stat
(SELECT SUM(sngl.Gamma2 * sngl.Gamma2) FROM sngl_inspiral as sngl WHERE sngl.event_id IN (SELECT event_id FROM coinc_event_map WHERE coinc_event_map.coinc_event_id == coinc_event.coinc_event_id)) AS ranking_stat
FROM
coinc_event
JOIN coinc_inspiral ON (
......
......@@ -182,7 +182,7 @@ CREATE TEMPORARY TABLE _cluster_info_ AS
(coinc_inspiral.end_time - (SELECT MIN(end_time) FROM coinc_inspiral)) + 1e-9 * coinc_inspiral.end_time_ns AS end_time,
coinc_inspiral.snr AS snr,
-- NOTE FIXME we are using bank_chisq to hold a chisq weighted snr value
(SELECT SUM(sngl.bank_chisq * sngl.bank_chisq) FROM sngl_inspiral as sngl WHERE sngl.event_id IN (SELECT event_id FROM coinc_event_map WHERE coinc_event_map.coinc_event_id == coinc_event.coinc_event_id)) AS ranking_stat
(SELECT SUM(sngl.Gamma2 * sngl.Gamma2) FROM sngl_inspiral as sngl WHERE sngl.event_id IN (SELECT event_id FROM coinc_event_map WHERE coinc_event_map.coinc_event_id == coinc_event.coinc_event_id)) AS ranking_stat
FROM
coinc_event
JOIN coinc_inspiral ON (
......
......@@ -47,6 +47,7 @@
#include <complex.h>
#include <math.h>
#include <string.h>
/*
......@@ -75,7 +76,7 @@
#include <gstlal_autocorrelation_chi2.h>
#include <gstlal_peakfinder.h>
/*
* ============================================================================
......@@ -476,3 +477,63 @@ unsigned gstlal_autocorrelation_chi2_float(
return output_length;
}
gsl_vector *gstlal_bankcorrelation_chi2_compute_norms(const gsl_matrix_complex *bankcorrelation_matrix)
{
gsl_vector *norm;
unsigned i;
unsigned channel;
unsigned channels = bankcorrelation_matrix->size1;
float complex cij;
float norms;
norm = gsl_vector_alloc(channels);
for(channel = 0; channel < channels; channel++) {
norms = 2*channels;
for(i = 0; i < channels; i++) {
cij = ((float) GSL_REAL(gsl_matrix_complex_get(bankcorrelation_matrix, i, channel)) + (float) GSL_IMAG(gsl_matrix_complex_get(bankcorrelation_matrix, i, channel)) * I);
norms -= 0.5*conjf(cij)*cij;
}
gsl_vector_set(norm, channel, norms);
}
return norm;
}
//FIXME make this like the float version
unsigned gstlal_bankcorrelation_chi2_from_peak(double *out, struct gstlal_peak_state *state, const gsl_matrix_complex *bcmat, const gsl_vector *bcnorm, const complex double *data, guint pad)
{
// FIXME add something
return 0;
}
unsigned gstlal_bankcorrelation_chi2_from_peak_float(float *out, struct gstlal_peak_state *state, const gsl_matrix_complex *bcmat, const gsl_vector *bcnorm, const complex float *data, guint pad)
{
unsigned i,j;
float complex *snr = state->values.as_float_complex;
float complex cij;
float complex xij;
float complex snrj;
int index;
for (i = 0; i < state->channels; i++)
{
out[i] = 0.0;
/*
* Don't bother computing if the event is below threshold, i.e.
* set to 0
*/
if (snr[i] == 0) continue;
for (j = 0; j < state->channels; j++)
{
index = (state->samples[i]) * state->channels + j;
snrj = *(data + index);
// Normalization of cij is in svd_bank.py. Phase must be (+).
cij = ((float) GSL_REAL(gsl_matrix_complex_get(bcmat, i, j)) + (float) GSL_IMAG(gsl_matrix_complex_get(bcmat, i, j)) * I);
// 0.5 is necessary due to the template normalization gstlal uses.
xij = 0.5 * snr[i] * cij - snrj;
out[i] += conjf(xij) * xij;
}
out[i] /= (float) gsl_vector_get(bcnorm, i);
}
return 0;
}
......@@ -44,6 +44,13 @@
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
/*
* Our own stuff
*/
#include <gstlal/gstlal_peakfinder.h>
/*
* ============================================================================
......@@ -59,3 +66,11 @@ unsigned gstlal_autocorrelation_chi2_autocorrelation_length(const gsl_matrix_com
gsl_vector *gstlal_autocorrelation_chi2_compute_norms(const gsl_matrix_complex *, const gsl_matrix_int *);
unsigned gstlal_autocorrelation_chi2(double *, const complex double *, unsigned, int, double, const gsl_matrix_complex *, const gsl_matrix_int *, const gsl_vector *);
unsigned gstlal_autocorrelation_chi2_float(float *, const float complex *, unsigned, int, double, const gsl_matrix_complex *, const gsl_matrix_int *, const gsl_vector *);
/*
* Bank veto functions
*/
gsl_vector *gstlal_bankcorrelation_chi2_compute_norms(const gsl_matrix_complex *bankcorrelation_matrix);
unsigned gstlal_bankcorrelation_chi2_from_peak(double *, struct gstlal_peak_state *, const gsl_matrix_complex *, const gsl_vector *, const complex double *, guint);
unsigned gstlal_bankcorrelation_chi2_from_peak_float(float *, struct gstlal_peak_state *, const gsl_matrix_complex *, const gsl_vector *, const complex float *, guint);
......@@ -141,6 +141,25 @@ int NAME(gstlal,series_around_peak)(struct gstlal_peak_state *state, TYPE *data,
return 0;
}
int NAME(gstlal,apply_threshold_to_peak)(struct gstlal_peak_state *state, double thresh)
{
guint channel;
TYPE *maxdata = MEMBER(state->values.as);
TYPE *interpmaxdata = MEMBER(state->interpvalues.as);
state->num_events = 0;
/* Decide if there are any events to keep */
for(channel = 0; channel < state->channels; channel++) {
if ( ABSFUNC(maxdata[channel]) <= thresh ) {
maxdata[channel] = 0;
interpmaxdata[channel] = 0;
state->samples[channel] = 0;
state->interpsamples[channel] = 0;
}
else state->num_events += 1;
}
return 0;
}
#undef PASTER3
#undef EVALUATOR3
#undef NAME
......
......@@ -14,6 +14,9 @@ int NAME(gstlal,series_around_peak)(struct gstlal_peak_state *state, TYPE *data,
/* fill the output */
int NAME(gstlal,fill_output_with_peak)(struct gstlal_peak_state *state, TYPE *data, guint64 length);
/* apply an additional threshold to an existing peak state */
int NAME(gstlal,apply_threshold_to_peak)(struct gstlal_peak_state *state, double thresh);
#undef PASTER
#undef EVALUATOR
#undef NAME
Markdown is supported
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