Loading [MathJax]/extensions/TeX/AMSmath.js
LALInspiral 5.0.3.1-ea7c608
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
_thinca.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2018 Kipp C. Cannon
3 *
4 * This program is free software; you can redistribute it and/or modify it
5 * under the terms of the GNU General Public License as published by the
6 * Free Software Foundation; either version 2 of the License, or (at your
7 * option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful, but
10 * WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 * General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License along
15 * with this program; if not, write to the Free Software Foundation, Inc.,
16 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17 */
18
19
20/*
21 * how does this work? this is a "get_coincs" class for the gstlal
22 * inspiral search to use with the snglcoinc machinery. a get_coincs
23 * object is a callable object initialized with a time-ordered sequence of
24 * triggers from one detector. subsequently, when called with a single
25 * trigger (from some other detector) it returns the in-order sequence of
26 * triggers from its initialization set that are coincident with the given
27 * trigger. this provides the double-coincidences upon which the rest of
28 * the coincidence machinery is built.
29 *
30 * for the gstlal inspiral search, two triggers are coincident if they were
31 * obtained from the same template and are within some time window of each
32 * other. to find coincidences quickly, the initial sequence of triggers
33 * is organized into a data structure like:
34 *
35 * array of event lists:
36 * template ID 0: [event, event, event, ...]
37 * template ID 1: [event, event, event, ...]
38 * ...
39 *
40 * the event lists are sorted by template ID, and the events in each list
41 * are time-ordered (which is ensured by providing them to the
42 * initialization code in time order). the events coincident with a
43 * trigger are found by retrieving the event list matching that trigger's
44 * ID with a bisection search, then the start and stop of the range of
45 * times coincident with the trigger are computed and the sequence of
46 * events that fall within that time range is returned by performing a
47 * bisection search to identify the first such event, and then a brute
48 * force search to find the last one (measured to be faster than a second
49 * bisection search).
50 */
51
52
53/*
54 * ============================================================================
55 *
56 * Preamble
57 *
58 * ============================================================================
59 */
60
61
62/* Ignore warnings in Python API itself */
63#pragma GCC diagnostic push
64#pragma GCC diagnostic ignored "-Wcast-align"
65#include <Python.h>
66#pragma GCC diagnostic pop
67
68#include <math.h>
69#include <stdlib.h>
70
71#include <lal/LALAtomicDatatypes.h>
72
73
74/*
75 * ============================================================================
76 *
77 * Internal Code
78 *
79 * ============================================================================
80 */
81
82
83/* retrieve the template ID of an event. in modern versions of the
84 * pipeline any integer may be used as a template ID, including negative
85 * values, so error status must be checked in the calling code using
86 * PyErr_Occurred() instead of by the return value alone. mirroring the
87 * behaviour of PyLong_AsLong(), we return -1 on error so that calling
88 * codes can test for that special value before calling PyErr_Occurred() to
89 * avoid the cost of that additional test when not needed */
90
91
92static long get_template_id(PyObject *event)
93{
94 PyObject *template_id = PyObject_GetAttrString(event, "template_id");
95 long val;
96
97 if(!template_id)
98 return -1;
99
100 val = PyLong_AsLong(template_id);
101 Py_DECREF(template_id);
102 /* PyLong_AsLong()'s error code is the same as ours, so we don't do
103 * error checking here. the calling code is responsible */
104
105 return val;
106}
107
108
109/* get the end time of an event as integer nanoseconds from the GPS epoch.
110 * returns -1 on error, but use PyErr_Occured() to confirm that an error
111 * has occured, and that -1 is not actually the return value. */
112
113
114static INT8 get_end_ns(PyObject *event)
115{
116 INT8 val = -1;
117 PyObject *end = PyObject_GetAttrString(event, "end");
118 if(end) {
119 PyObject *ns = PyObject_CallMethod(end, "ns", NULL);
120 Py_DECREF(end);
121 if(ns) {
122 val = PyLong_AsLong(ns);
123 Py_DECREF(ns);
124 }
125 }
126 return val;
127}
128
129
130/* compute the range of times, [start, stop), for which triggers from other
131 * detectors are coincident with this one. */
132
133
134static int compute_start_stop(INT8 *start_ns, INT8 *stop_ns, PyObject *event, PyObject *offset, PyObject* coinc_window)
135{
136 INT8 end_ns = get_end_ns(event);
137 double offset_ns = PyFloat_AsDouble(offset);
138 double coinc_window_ns = PyFloat_AsDouble(coinc_window);
139
140 if((end_ns == -1 || offset_ns == -1. || coinc_window_ns == -1.) && PyErr_Occurred())
141 return -1;
142 offset_ns *= 1e9;
143 coinc_window_ns *= 1e9;
144
145 *start_ns = end_ns + lround(offset_ns - coinc_window_ns);
146 *stop_ns = end_ns + lround(offset_ns + coinc_window_ns);
147
148 return 0;
149}
150
151
152/*
153 * event table
154 */
155
156
157struct event {
159 PyObject *event;
160};
161
162
165 struct event *events;
167};
168
169/* comparison operator used with bsearch() to find an event list
170 * corresponding to a given template ID */
171
172static int event_sequence_get_cmp(const void *key, const void *obj)
173{
174 return (long) key - ((const struct event_sequence *) obj)->template_id;
175}
176
177/* comparison operator used with qsort() to sort event lists by template
178 * ID */
179
180static int event_sequence_sort_cmp(const void *a, const void *b)
181{
182 return ((const struct event_sequence *) a)->template_id - ((const struct event_sequence *) b)->template_id;
183}
184
185/* bsearch() wrapper to retrieve an event list by template ID */
186
187static struct event_sequence *event_sequence_get(struct event_sequence *event_sequences, int n, long template_id)
188{
189 return bsearch((void *) template_id, event_sequences, n, sizeof(*event_sequences), event_sequence_get_cmp);
190}
191
192/* insert a time-ordered event into the array of event list arrays */
193
194static int event_sequence_insert(struct event_sequence **event_sequences, int *n, PyObject *event)
195{
197 INT8 end_ns = get_end_ns(event);
199 int need_sort = 0;
200
201 if((template_id == -1 || end_ns == -1) && PyErr_Occurred())
202 return -1;
203
204 /* get the event list for this template ID or create a new empty
205 * one if there isn't one yet. the new list is appeneded, and the
206 * array of event lists put into the correct order later. */
207
208 event_sequence = event_sequence_get(*event_sequences, *n, template_id);
209 if(!event_sequence) {
210 struct event_sequence *new = realloc(*event_sequences, (*n + 1) * sizeof(**event_sequences));
211 if(!new) {
212 PyErr_SetString(PyExc_MemoryError, "realloc() failed");
213 return -1;
214 }
215 *event_sequences = new;
216 event_sequence = &new[*n];
217 (*n)++;
219 event_sequence->events = NULL;
221 need_sort = 1;
222 }
223
224 /* make room for the new event. allocate space for triggers 128 at
225 * a time */
226
227 if(!(event_sequence->length & 0x7f)) {
228 struct event *new = realloc(event_sequence->events, ((event_sequence->length | 0x7f) + 1) * sizeof(*event_sequence->events));
229 if(!new) {
230 PyErr_SetString(PyExc_MemoryError, "realloc() failed");
231 return -1;
232 }
233 event_sequence->events = new;
234 }
235
236 /* record the event in the event list */
237
238 Py_INCREF(event);
242
243 /* if a new event list was created, put the array of event lists in
244 * order by template ID */
245
246 if(need_sort)
247 qsort(*event_sequences, *n, sizeof(**event_sequences), event_sequence_sort_cmp);
248
249 return 0;
250}
251
252
253static void event_sequence_free(struct event_sequence *event_sequences, int n)
254{
255 while(n--) {
256 struct event_sequence *event_sequence = &event_sequences[n];
257 while(event_sequence->length--)
259 free(event_sequence->events);
260 event_sequence->events = NULL;
261 }
262
263 free(event_sequences);
264}
265
266
267static Py_ssize_t bisect_left(struct event *items, Py_ssize_t hi, INT8 end_ns)
268{
269 Py_ssize_t lo = 0;
270
271 while(lo < hi) {
272 Py_ssize_t mid = (lo + hi) / 2;
273 if(items[mid].end_ns < end_ns)
274 lo = mid + 1;
275 else
276 /* items[mid].end_ns >= end_ns */
277 hi = mid;
278 }
279
280 return lo;
281}
282
283static PyObject *event_sequence_extract(struct event_sequence *event_sequence, INT8 start_ns, INT8 stop_ns)
284{
285 Py_ssize_t lo = bisect_left(event_sequence->events, event_sequence->length, start_ns);
286 Py_ssize_t hi;
287 Py_ssize_t i;
288 PyObject *result;
289
290 if(lo < 0)
291 return NULL;
292
293 for(hi = lo; hi < event_sequence->length && event_sequence->events[hi].end_ns < stop_ns; hi++);
294
295 result = PyTuple_New(hi - lo);
296 if(!result)
297 return NULL;
298 for(i = 0; lo < hi; i++, lo++) {
299 PyObject *item = event_sequence->events[lo].event;
300 Py_INCREF(item);
301 PyTuple_SET_ITEM(result, i, item);
302 }
303
304 return result;
305}
306
307
308/*
309 * ============================================================================
310 *
311 * get_coincs Class
312 *
313 * ============================================================================
314 */
315
316
317/*
318 * Structure
319 */
320
321
323 PyObject_HEAD
324
327};
328
329
330/*
331 * __del__() method
332 */
333
334
335static void get_coincs__del__(PyObject *self)
336{
337 struct get_coincs *get_coincs = (struct get_coincs *) self;
338
342
343 self->ob_type->tp_free(self);
344}
345
346
347/*
348 * __init__() method
349 */
350
351
352static int get_coincs__init__(PyObject *self, PyObject *args, PyObject *kwds)
353{
354 struct get_coincs *get_coincs = (struct get_coincs *) self;
355 PyObject *events;
356 PyObject **item, **last;
357
358 if(kwds) {
359 PyErr_SetString(PyExc_ValueError, "unexpected keyword arguments");
360 return -1;
361 }
362 if(!PyArg_ParseTuple(args, "O", &events))
363 return -1;
364
369
370 item = PySequence_Fast_ITEMS(events);
371 if(!item)
372 return -1;
373 last = item + PySequence_Length(events);
374 for(; item < last; item++)
379 return -1;
380 }
381
382 return 0;
383}
384
385
386/*
387 * __call__() method
388 */
389
390
391static PyObject *get_coincs__call__(PyObject *self, PyObject *args, PyObject *kwds)
392{
393 struct get_coincs *get_coincs = (struct get_coincs *) self;
394 PyObject *event_a;
395 PyObject *offset_a;
396 PyObject *coinc_window;
397 long template_id;
399 INT8 start_ns, stop_ns;
400
401 if(kwds) {
402 PyErr_SetString(PyExc_ValueError, "unexpected keyword arguments");
403 return NULL;
404 }
405 if(!PyArg_ParseTuple(args, "OOO", &event_a, &offset_a, &coinc_window))
406 return NULL;
407
408 template_id = get_template_id(event_a);
409 if(template_id == -1 && PyErr_Occurred())
410 return NULL;
413 /* no events */
414 Py_INCREF(Py_None);
415 return Py_None;
416 }
417
418 if(compute_start_stop(&start_ns, &stop_ns, event_a, offset_a, coinc_window) < 0)
419 return NULL;
420
421 /* return the events matching event_a's .template_id and having end
422 * times in the range t - coinc_window <= end < t + coinc_window
423 * where t = event_a.end + offset_a */
424
425 return event_sequence_extract(event_sequence, start_ns, stop_ns);
426}
427
428
429/*
430 * Type information
431 */
432
433
434static PyTypeObject get_coincs_Type = {
435 PyVarObject_HEAD_INIT(NULL, 0)
436 .tp_name = MODULE_NAME ".get_coincs",
437 .tp_basicsize = sizeof(struct get_coincs),
438 .tp_dealloc = get_coincs__del__,
439 .tp_call = get_coincs__call__,
440 .tp_doc = "",
441 .tp_flags = Py_TPFLAGS_DEFAULT,
442 .tp_init = get_coincs__init__,
443 .tp_new = PyType_GenericNew,
444};
445
446
447/*
448 * ============================================================================
449 *
450 * Entry Point
451 *
452 * ============================================================================
453 */
454
455
456PyMODINIT_FUNC PyInit__thinca(void); /* silence warning */
457PyMODINIT_FUNC PyInit__thinca(void)
458{
459 static struct PyModuleDef modef = {
460 PyModuleDef_HEAD_INIT,
461 .m_name = MODULE_NAME,
462 .m_doc = "",
463 .m_size = -1,
464 .m_methods = NULL,
465 };
466 PyObject *module = PyModule_Create(&modef);
467
468 if(PyType_Ready(&get_coincs_Type) < 0) {
469 Py_DECREF(module);
470 module = NULL;
471 goto done;
472 }
473 Py_INCREF((PyObject *) &get_coincs_Type);
474 if(PyModule_AddObject(module, "get_coincs", (PyObject *) &get_coincs_Type) < 0) {
475 Py_DECREF(module);
476 module = NULL;
477 goto done;
478 }
479
480done:
481 return module;
482}
static Py_ssize_t bisect_left(struct event *items, Py_ssize_t hi, INT8 end_ns)
Definition: _thinca.c:267
static PyObject * event_sequence_extract(struct event_sequence *event_sequence, INT8 start_ns, INT8 stop_ns)
Definition: _thinca.c:283
static int event_sequence_insert(struct event_sequence **event_sequences, int *n, PyObject *event)
Definition: _thinca.c:194
static INT8 get_end_ns(PyObject *event)
Definition: _thinca.c:114
static int event_sequence_get_cmp(const void *key, const void *obj)
Definition: _thinca.c:172
static void event_sequence_free(struct event_sequence *event_sequences, int n)
Definition: _thinca.c:253
PyMODINIT_FUNC PyInit__thinca(void)
Definition: _thinca.c:457
static PyObject * get_coincs__call__(PyObject *self, PyObject *args, PyObject *kwds)
Definition: _thinca.c:391
static long get_template_id(PyObject *event)
Definition: _thinca.c:92
static struct event_sequence * event_sequence_get(struct event_sequence *event_sequences, int n, long template_id)
Definition: _thinca.c:187
static int event_sequence_sort_cmp(const void *a, const void *b)
Definition: _thinca.c:180
static void get_coincs__del__(PyObject *self)
Definition: _thinca.c:335
static int get_coincs__init__(PyObject *self, PyObject *args, PyObject *kwds)
Definition: _thinca.c:352
static PyTypeObject get_coincs_Type
Definition: _thinca.c:434
static int compute_start_stop(INT8 *start_ns, INT8 *stop_ns, PyObject *event, PyObject *offset, PyObject *coinc_window)
Definition: _thinca.c:134
double i
int64_t INT8
static const INT4 a
int n
struct event * events
Definition: _thinca.c:165
long template_id
Definition: _thinca.c:164
Definition: _thinca.c:157
INT8 end_ns
Definition: _thinca.c:158
PyObject * event
Definition: _thinca.c:159
PyObject_HEAD struct event_sequence * event_sequences
Definition: _thinca.c:325
int n_sequences
Definition: _thinca.c:326