Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
ResultsToplist.c
Go to the documentation of this file.
1//
2// Copyright (C) 2016, 2017 Karl Wette
3//
4// This program is free software; you can redistribute it and/or modify
5// it under the terms of the GNU General Public License as published by
6// the Free Software Foundation; either version 2 of the License, or
7// (at your option) any later version.
8//
9// This program is distributed in the hope that it will be useful,
10// but WITHOUT ANY WARRANTY; without even the implied warranty of
11// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12// GNU General Public License for more details.
13//
14// You should have received a copy of the GNU General Public License
15// along with with program; see the file COPYING. If not, write to the
16// Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17// MA 02110-1301 USA
18//
19
20///
21/// \file
22/// \ingroup lalpulsar_bin_Weave
23///
24
25#include "ResultsToplist.h"
26#include "ComputeResults.h"
27
28#include <lal/VectorMath.h>
29#include <lal/UserInputPrint.h>
30
31// Compare two quantities, and return a sort order value if they are unequal
32#define COMPARE_BY( x, y ) do { if ( (x) < (y) ) return -1; if ( (x) > (y) ) return +1; } while(0)
33
34///
35/// Toplist of output results
36///
38 /// Struct holding all parameters for which statistics to output and compute, when, and how
39 /// NOTE: this is only a reference-pointer to WeaveStatisticsParams, while WeaveSemiResults is the *owner*
40 WeaveStatisticsParams *statistics_params;
41 /// Number of spindown parameters to output
42 size_t nspins;
43 /// Name of ranking statistic
44 const char *stat_name;
45 /// Description of ranking statistic
46 const char *stat_desc;
47 /// Function which returns pointer to array of statistics by which toplist items are ranked
49 /// Function which returns the value of the statistic by which toplist items are ranked
51 /// Function which sets the value of the statistic by which toplist items are ranked
53 /// Vector of indexes of toplist results which should be considered for addition
55 /// Serial number which is incremented for each item
57 /// Heap which ranks toplist items
58 LALHeap *heap;
59 /// Save a no-longer-used toplist item for re-use
60 WeaveResultsToplistItem *saved_item;
61};
62
63///
64/// \name Internal functions
65///
66/// @{
67
68static WeaveResultsToplistItem *toplist_item_create( const WeaveResultsToplist *toplist );
69static int compare_templates( BOOLEAN *equal, const char *loc_str, const char *tmpl_str, const UINT4 round_param_to_dp, const UINT4 round_param_to_sf, const REAL8 param_tol_mism, const gsl_matrix *metric, const SuperskyTransformData *rssky_transf, const UINT8 serial_1, const UINT8 serial_2, const PulsarDopplerParams *phys_1, const PulsarDopplerParams *phys_2 );
70static int compare_vectors( BOOLEAN *equal, const VectorComparison *result_tol, const REAL4Vector *res_1, const REAL4Vector *res_2, const UINT4 r1, const UINT4 r2 );
71static double round_to_dp_sf( double x, const UINT4 dp, const UINT4 sf );
72static int toplist_fits_table_init( FITSFile *file, const WeaveResultsToplist *toplist );
73static int toplist_fits_table_write_visitor( void *param, const void *x );
74static int toplist_item_sort_by_semi_phys( const void *x, const void *y );
75static int toplist_item_sort_by_serial( const void *x, const void *y );
76static void toplist_item_destroy( WeaveResultsToplistItem *item );
77static int toplist_item_compare( void *param, const void *x, const void *y );
78static int toplist_fill_completionloop_stats( void *param, void *x );
79
80/// @}
81
82///
83/// Round a floating-point number 'x' to 'dp' decimal places (if > 0),
84/// then 'sf' significant figures (if > 0).
85///
86static double round_to_dp_sf(
87 double x,
88 const UINT4 dp,
89 const UINT4 sf
90)
91{
92 char buf[32];
93 if ( dp > 0 ) {
94 snprintf( buf, sizeof( buf ), "%0.*f", dp, x );
95 x = atof( buf );
96 }
97 if ( sf > 0 ) {
98 snprintf( buf, sizeof( buf ), "%0.*e", sf - 1, x );
99 x = atof( buf );
100 }
101 return x;
102}
103
104///
105/// Create a toplist item
106///
107WeaveResultsToplistItem *toplist_item_create(
108 const WeaveResultsToplist *toplist
109)
110{
111
112 // Check input
113 XLAL_CHECK_NULL( toplist != NULL, XLAL_EFAULT );
114
115 // Allocate memory for item
116 WeaveResultsToplistItem *item = XLALCalloc( 1, sizeof( *item ) );
117 XLAL_CHECK_NULL( item != NULL, XLAL_ENOMEM );
118
119 const WeaveStatisticsParams *params = toplist-> statistics_params;
120 WeaveStatisticType store_per_segment_stats = ( params->statistics_to_output[0] & ( WEAVE_STATISTIC_COH2F | WEAVE_STATISTIC_COH2F_DET ) );
121 WeaveStatisticType store_coh2F[2];
122 WeaveStatisticType store_coh2F_det[2];
123 store_coh2F[0] = ( params->mainloop_statistics_to_keep & WEAVE_STATISTIC_COH2F );
124 store_coh2F_det[0] = ( params->mainloop_statistics_to_keep & WEAVE_STATISTIC_COH2F_DET );
125 store_coh2F[1] = ( params->completionloop_statistics[1] & WEAVE_STATISTIC_COH2F );
126 store_coh2F_det[1] = ( params->completionloop_statistics[1] & WEAVE_STATISTIC_COH2F_DET );
127
128 // Allocate memory for per-segment output results
129 if ( store_per_segment_stats ) {
130 item->coh_alpha = XLALCalloc( params->nsegments, sizeof( *item->coh_alpha ) );
131 XLAL_CHECK_NULL( item->coh_alpha != NULL, XLAL_ENOMEM );
132 item->coh_delta = XLALCalloc( params->nsegments, sizeof( *item->coh_delta ) );
133 XLAL_CHECK_NULL( item->coh_delta != NULL, XLAL_ENOMEM );
134 for ( size_t k = 0; k <= toplist->nspins; ++k ) {
135 item->coh_fkdot[k] = XLALCalloc( params->nsegments, sizeof( *item->coh_fkdot[k] ) );
136 XLAL_CHECK_NULL( item->coh_fkdot[k] != NULL, XLAL_ENOMEM );
137 }
138 }
139
140 // Allocate memory for per-segment coh2F and coh2F_det statistics if requested
141 for ( UINT4 istage = 0; istage < 2; istage ++ ) {
142 if ( store_coh2F[istage] ) {
143 item->stage[istage].coh2F = XLALCalloc( params->nsegments, sizeof( *item->stage[istage].coh2F ) );
144 XLAL_CHECK_NULL( item->stage[istage].coh2F != NULL, XLAL_ENOMEM );
145 }
146 if ( store_coh2F_det[istage] ) {
147 for ( UINT4 X = 0; X < params->detectors->length; ++X ) {
148 item->stage[istage].coh2F_det[X] = XLALCalloc( params->nsegments, sizeof( *item->stage[istage].coh2F_det[X] ) );
149 XLAL_CHECK_NULL( item->stage[istage].coh2F_det[X] != NULL, XLAL_ENOMEM );
150 }
151 }
152 }
153
154 return item;
155
156}
157
158///
159/// Destroy a toplist item
160///
162 WeaveResultsToplistItem *item
163)
164{
165 if ( item != NULL ) {
166 XLALFree( item->coh_alpha );
167 XLALFree( item->coh_delta );
168 for ( size_t k = 0; k < PULSAR_MAX_SPINS; ++k ) {
169 XLALFree( item->coh_fkdot[k] );
170 }
171 for ( size_t istage = 0; istage < 2; ++istage ) {
172 XLALFree( item->stage[istage].coh2F );
173 for ( size_t X = 0; X < PULSAR_MAX_DETECTORS; ++X ) {
174 XLALFree( item->stage[istage].coh2F_det[X] );
175 }
176 }
177 XLALFree( item );
178 }
179}
180
181///
182/// Compare toplist items
183///
185 void *param,
186 const void *x,
187 const void *y
188)
189{
191 const WeaveResultsToplistItem *ix = ( const WeaveResultsToplistItem * ) x;
192 const WeaveResultsToplistItem *iy = ( const WeaveResultsToplistItem * ) y;
193 COMPARE_BY( item_get_rank_stat_fcn( iy ), item_get_rank_stat_fcn( ix ) ); // Compare in descending order
194 return 0;
195}
196
197///
198/// Initialise a FITS table for writing/reading a toplist
199///
201 FITSFile *file,
202 const WeaveResultsToplist *toplist
203)
204{
205
206 // Check input
207 XLAL_CHECK( file != NULL, XLAL_EFAULT );
208
209 char col_name[64];
210
211 // Begin FITS table description
212 XLAL_FITS_TABLE_COLUMN_BEGIN( WeaveResultsToplistItem );
213
214 // Add columns for semicoherent template parameters
216 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_NAMED( file, REAL8, semi_alpha, "alpha [rad]" ) == XLAL_SUCCESS, XLAL_EFUNC );
217 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_NAMED( file, REAL8, semi_delta, "delta [rad]" ) == XLAL_SUCCESS, XLAL_EFUNC );
218 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_NAMED( file, REAL8, semi_fkdot[0], "freq [Hz]" ) == XLAL_SUCCESS, XLAL_EFUNC );
219 for ( size_t k = 1; k <= toplist->nspins; ++k ) {
220 snprintf( col_name, sizeof( col_name ), "f%zudot [Hz/s^%zu]", k, k );
222 }
223
224 const WeaveStatisticsParams *params = toplist->statistics_params;
225
226 const char *stage_suffix[2] = {"", "_rec"};
227
228 for ( UINT4 istage = 0; istage < 2; istage ++ ) {
229 // Which statistics have been requested for output at this stage ('stage0' vs 'recalc')?
230 WeaveStatisticType statistics_to_output = params->statistics_to_output[istage];
231
232 // Add column for mean multi-detector F-statistic
233 if ( statistics_to_output & WEAVE_STATISTIC_MEAN2F ) {
234 snprintf( col_name, sizeof( col_name ), "%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_MEAN2F ), stage_suffix[istage] );
236 }
237
238 // Add columns for mean per-detector F-statistic
239 if ( statistics_to_output & WEAVE_STATISTIC_MEAN2F_DET ) {
240 for ( size_t i = 0; i < params->detectors->length; ++i ) {
241 snprintf( col_name, sizeof( col_name ), "%s_%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_MEAN2F ), params->detectors->data[i], stage_suffix[istage] );
242 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_NAMED( file, REAL4, stage[istage].mean2F_det[i], col_name ) == XLAL_SUCCESS, XLAL_EFUNC );
243 }
244 }
245
246 // Add column for multi-detector max2F statistic
247 if ( statistics_to_output & WEAVE_STATISTIC_MAX2F ) {
248 snprintf( col_name, sizeof( col_name ), "%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_MAX2F ), stage_suffix[istage] );
250 }
251
252 // Add columns for per-detector max2F_det statistic
253 if ( statistics_to_output & WEAVE_STATISTIC_MAX2F_DET ) {
254 for ( size_t i = 0; i < params->detectors->length; ++i ) {
255 snprintf( col_name, sizeof( col_name ), "%s_%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_MAX2F ), params->detectors->data[i], stage_suffix[istage] );
256 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_NAMED( file, REAL4, stage[istage].max2F_det[i], col_name ) == XLAL_SUCCESS, XLAL_EFUNC );
257 }
258 }
259
260 // Add column for multi-detector sum2F statistic
261 if ( statistics_to_output & WEAVE_STATISTIC_SUM2F ) {
262 snprintf( col_name, sizeof( col_name ), "%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_SUM2F ), stage_suffix[istage] );
264 }
265
266 // Add columns for per-detector sum2F statistic
267 if ( statistics_to_output & WEAVE_STATISTIC_SUM2F_DET ) {
268 for ( size_t i = 0; i < params->detectors->length; ++i ) {
269 snprintf( col_name, sizeof( col_name ), "%s_%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_SUM2F ), params->detectors->data[i], stage_suffix[istage] );
270 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_NAMED( file, REAL4, stage[istage].sum2F_det[i], col_name ) == XLAL_SUCCESS, XLAL_EFUNC );
271 }
272 }
273
274 // Add column for BSGL statistic
275 if ( statistics_to_output & WEAVE_STATISTIC_BSGL ) {
276 snprintf( col_name, sizeof( col_name ), "%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_BSGL ), stage_suffix[istage] );
277 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_NAMED( file, REAL4, stage[istage].log10BSGL, col_name ) == XLAL_SUCCESS, XLAL_EFUNC );
278 }
279
280 // Add column for BSGLtL statistic
281 if ( statistics_to_output & WEAVE_STATISTIC_BSGLtL ) {
282 snprintf( col_name, sizeof( col_name ), "%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_BSGLtL ), stage_suffix[istage] );
283 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_NAMED( file, REAL4, stage[istage].log10BSGLtL, col_name ) == XLAL_SUCCESS, XLAL_EFUNC );
284 }
285
286 // Add column for BtSGLtL statistic
287 if ( statistics_to_output & WEAVE_STATISTIC_BtSGLtL ) {
288 snprintf( col_name, sizeof( col_name ), "%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_BtSGLtL ), stage_suffix[istage] );
289 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_NAMED( file, REAL4, stage[istage].log10BtSGLtL, col_name ) == XLAL_SUCCESS, XLAL_EFUNC );
290 }
291
292 // Add column for multi-detector number-count statistic
293 if ( statistics_to_output & WEAVE_STATISTIC_NCOUNT ) {
294 snprintf( col_name, sizeof( col_name ), "%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_NCOUNT ), stage_suffix[istage] );
296 }
297
298 // Add column for per-detector number-count statistics
299 if ( statistics_to_output & WEAVE_STATISTIC_NCOUNT_DET ) {
300 for ( size_t i = 0; i < params->detectors->length; ++i ) {
301 snprintf( col_name, sizeof( col_name ), "%s_%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_NCOUNT ), params->detectors->data[i], stage_suffix[istage] );
302 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_NAMED( file, REAL4, stage[istage].ncount_det[i], col_name ) == XLAL_SUCCESS, XLAL_EFUNC );
303 }
304 }
305
306 //
307 // Add columns for per-segment statistics
308 //
309
310 if ( istage == 0 ) { // Only output these for 'stage 0'
311 // We tie the output of per-segment coordinates to the output of any per-segment statistics
312 WeaveStatisticType per_segment_stats = ( statistics_to_output & ( WEAVE_STATISTIC_COH2F | WEAVE_STATISTIC_COH2F_DET ) );
313 // Add columns for coherent template parameters
314 if ( per_segment_stats ) {
315 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_PTR_ARRAY_NAMED( file, REAL8, params->nsegments, coh_alpha, "alpha_seg [rad]" ) == XLAL_SUCCESS, XLAL_EFUNC );
316 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_PTR_ARRAY_NAMED( file, REAL8, params->nsegments, coh_delta, "delta_seg [rad]" ) == XLAL_SUCCESS, XLAL_EFUNC );
317 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_PTR_ARRAY_NAMED( file, REAL8, params->nsegments, coh_fkdot[0], "freq_seg [Hz]" ) == XLAL_SUCCESS, XLAL_EFUNC );
318 for ( size_t k = 1; k <= toplist->nspins; ++k ) {
319 snprintf( col_name, sizeof( col_name ), "f%zudot_seg [Hz/s^%zu]", k, k );
321 }
322 }
323 }
324
325 // Add column for per-segment coherent multi-detector 2F statistics
326 if ( statistics_to_output & WEAVE_STATISTIC_COH2F ) {
327 snprintf( col_name, sizeof( col_name ), "%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_COH2F ), stage_suffix[istage] );
329 }
330
331 // Add columns for per-segment coherent per-detector 2F statistics
332 if ( statistics_to_output & WEAVE_STATISTIC_COH2F_DET ) {
333 for ( size_t i = 0; i < params->detectors->length; ++i ) {
334 snprintf( col_name, sizeof( col_name ), "%s_%s%s", WEAVE_STATISTIC_NAME( WEAVE_STATISTIC_COH2F ), params->detectors->data[i], stage_suffix[istage] );
335 XLAL_CHECK( XLAL_FITS_TABLE_COLUMN_ADD_PTR_ARRAY_NAMED( file, REAL4, params->nsegments, stage[istage].coh2F_det[i], col_name ) == XLAL_SUCCESS, XLAL_EFUNC );
336 }
337 }
338
339 }
340
341 return XLAL_SUCCESS;
342
343}
344
345///
346/// Function to update given toplist item with missing 'completion loop' statistics
347///
349 void *param,
350 void *x
351)
352{
353 XLAL_CHECK( param != NULL, XLAL_EFAULT );
354 XLAL_CHECK( x != NULL, XLAL_EFAULT );
355
356 WeaveResultsToplistItem *item = ( WeaveResultsToplistItem * ) x;
357 WeaveStatisticsParams *stats_params = ( WeaveStatisticsParams * )param;
358 UINT4 ndetectors = stats_params->detectors->length;
359 UINT4 nsegments = stats_params->nsegments;
360
361 for ( UINT4 istage = 0; istage < 2; istage ++ ) {
362 WeaveStatisticType stage_stats = stats_params->completionloop_statistics[istage];
363
364 // Re-calculate per-sement coherent 2F|2F_det statistics in semi-coherent template point
365 if ( stage_stats & ( WEAVE_STATISTIC_COH2F | WEAVE_STATISTIC_COH2F_DET ) ) {
366 XLAL_CHECK( istage > 0, XLAL_EERR, "BUG: requested 'coh2F' or 'coh2F_det' in stage0 completion loop ==> should never happen!\n" );
368 semi_phys.refTime = stats_params->ref_time;
369 semi_phys.Alpha = item->semi_alpha;
370 semi_phys.Delta = item->semi_delta;
371 memcpy( semi_phys.fkdot, item->semi_fkdot, sizeof( semi_phys.fkdot ) );
372
373 const UINT4 nfreqs = 1;
374 for ( size_t l = 0; l < nsegments; ++l ) {
375 XLAL_CHECK( XLALWeaveCohResultsCompute( &( stats_params->coh_res ), stats_params->coh_input_recalc[l], &semi_phys, nfreqs, NULL ) == XLAL_SUCCESS, XLAL_EFUNC );
376 REAL4Vector *coh2F = NULL;
378 BOOLEAN have_coh2F_det;
379 XLAL_CHECK( XLALWeaveCohResultsExtract( &coh2F, coh2F_det, &have_coh2F_det, stats_params->coh_res, stats_params->coh_input_recalc[l] ) == XLAL_SUCCESS, XLAL_EFUNC );
380 if ( stage_stats & WEAVE_STATISTIC_COH2F ) {
381 item->stage[istage].coh2F[l] = coh2F->data[0];
382 }
383 if ( stage_stats & WEAVE_STATISTIC_COH2F_DET ) {
384 XLAL_CHECK( have_coh2F_det, XLAL_EFAILED, "BUG: caller requested per-detector 2F-statistic, but was not computed\n" );
385 for ( UINT4 X = 0; X < ndetectors; ++X ) {
386 item->stage[istage].coh2F_det[X][l] = ( coh2F_det[X] != NULL ) ? coh2F_det[X]->data[0] : NAN;
387 }
388 }
389 }
390
391 }
392
393 if ( stage_stats & WEAVE_STATISTIC_MAX2F ) {
394 item->stage[istage].max2F = 0;
395 for ( size_t l = 0; l < nsegments; ++l ) {
396 item->stage[istage].max2F = fmaxf( item->stage[istage].max2F, item->stage[istage].coh2F[l] );
397 }
398 }
399
400 if ( stage_stats & WEAVE_STATISTIC_MAX2F_DET ) {
401 for ( size_t X = 0; X < ndetectors; ++X ) {
402 item->stage[istage].max2F_det[X] = 0;
403 for ( size_t l = 0; l < nsegments; ++l ) {
404 REAL4 item_Xl = item->stage[istage].coh2F_det[X][l];
405 item->stage[istage].max2F_det[X] = isnan( item_Xl ) ? item->stage[istage].max2F_det[X] : fmaxf( item->stage[istage].max2F_det[X], item_Xl );
406 }
407 }
408 }
409
410 if ( stage_stats & WEAVE_STATISTIC_SUM2F ) {
411 item->stage[istage].sum2F = 0;
412 for ( size_t l = 0; l < nsegments; ++l ) {
413 item->stage[istage].sum2F += item->stage[istage].coh2F[l];
414 }
415 }
416
417 if ( stage_stats & WEAVE_STATISTIC_SUM2F_DET ) {
418 for ( size_t X = 0; X < ndetectors; ++X ) {
419 item->stage[istage].sum2F_det[X] = 0;
420 for ( size_t l = 0; l < nsegments; ++l ) {
421 REAL4 item_Xl = item->stage[istage].coh2F_det[X][l];
422 item->stage[istage].sum2F_det[X] += isnan( item_Xl ) ? 0 : item_Xl;
423 }
424 }
425 }
426
427 if ( stage_stats & WEAVE_STATISTIC_MEAN2F ) {
428 item->stage[istage].mean2F = item->stage[istage].sum2F / nsegments;
429 }
430
431 if ( stage_stats & WEAVE_STATISTIC_MEAN2F_DET ) {
432 for ( size_t X = 0; X < ndetectors; ++X ) {
433 item->stage[istage].mean2F_det[X] = item->stage[istage].sum2F_det[X] / stats_params->n2F_det[X];
434 }
435 }
436
437 if ( stage_stats & WEAVE_STATISTIC_BSGL ) {
438 item->stage[istage].log10BSGL = XLALComputeBSGL( item->stage[istage].sum2F, item->stage[istage].sum2F_det, stats_params->BSGL_setup );
439 }
440
441 if ( stage_stats & WEAVE_STATISTIC_BSGLtL ) {
442 item->stage[istage].log10BSGLtL = XLALComputeBSGLtL( item->stage[istage].sum2F, item->stage[istage].sum2F_det, item->stage[istage].max2F_det, stats_params->BSGL_setup );
443 }
444
445 if ( stage_stats & WEAVE_STATISTIC_BtSGLtL ) {
446 item->stage[istage].log10BtSGLtL = XLALComputeBtSGLtL( item->stage[istage].max2F, item->stage[istage].sum2F_det, item->stage[istage].max2F_det, stats_params->BSGL_setup );
447 }
448
449 if ( stage_stats & WEAVE_STATISTIC_NCOUNT ) {
450 item->stage[istage].ncount = 0;
451 for ( size_t l = 0; l < nsegments; ++l ) {
452 item->stage[istage].ncount += ( item->stage[istage].coh2F[l] > stats_params->nc_2Fth ) ? 1 : 0;
453 }
454 }
455
456 if ( stage_stats & WEAVE_STATISTIC_NCOUNT_DET ) {
457 for ( size_t X = 0; X < ndetectors; ++X ) {
458 item->stage[istage].ncount_det[X] = 0;
459 for ( size_t l = 0; l < nsegments; ++l ) {
460 REAL4 item_Xl = item->stage[istage].coh2F_det[X][l];
461 if ( isnan( item_Xl ) ) {
462 continue;
463 }
464 item->stage[istage].ncount_det[X] += ( item_Xl > stats_params->nc_2Fth ) ? 1 : 0;
465 }
466 }
467 }
468
469 }
470
471 return XLAL_SUCCESS;
472
473}
474
475///
476/// Visitor function for writing a toplist to a FITS table
477///
479 void *param,
480 const void *x
481)
482{
483 FITSFile *file = ( FITSFile * ) param;
485 return XLAL_SUCCESS;
486}
487
488///
489/// Sort toplist items by physical coordinates of semicoherent template
490///
491/// For stable comparisons, the order of parameter comparisons should be the same
492/// as the order in which parameters are generated by the search lattice tiling.
493///
495 const void *x,
496 const void *y
497)
498{
499 const WeaveResultsToplistItem *ix = *( const WeaveResultsToplistItem * const * ) x;
500 const WeaveResultsToplistItem *iy = *( const WeaveResultsToplistItem * const * ) y;
501 COMPARE_BY( ix->semi_alpha, iy->semi_alpha ); // Compare in ascending order
502 COMPARE_BY( ix->semi_delta, iy->semi_delta ); // Compare in ascending order
503 for ( size_t s = 1; s < XLAL_NUM_ELEM( ix->semi_fkdot ); ++s ) {
504 COMPARE_BY( ix->semi_fkdot[s], iy->semi_fkdot[s] ); // Compare in ascending order
505 }
506 COMPARE_BY( ix->semi_fkdot[0], iy->semi_fkdot[0] ); // Compare in ascending order
507 return 0;
508}
509
510///
511/// Sort toplist items by serial number of template
512///
514 const void *x,
515 const void *y
516)
517{
518 const WeaveResultsToplistItem *ix = *( const WeaveResultsToplistItem * const * ) x;
519 const WeaveResultsToplistItem *iy = *( const WeaveResultsToplistItem * const * ) y;
520 COMPARE_BY( ix->serial, iy->serial ); // Compare in ascending order
521 return 0;
522}
523
524///
525/// Compute two template parameters
526///
528 BOOLEAN *equal,
529 const char *loc_str,
530 const char *tmpl_str,
531 const UINT4 round_param_to_dp,
532 const UINT4 round_param_to_sf,
533 const REAL8 param_tol_mism,
534 const gsl_matrix *metric,
535 const SuperskyTransformData *rssky_transf,
536 const UINT8 serial_1,
537 const UINT8 serial_2,
538 const PulsarDopplerParams *phys_1,
539 const PulsarDopplerParams *phys_2
540)
541{
542
543 // Check input
544 XLAL_CHECK( equal != NULL, XLAL_EINVAL );
545 XLAL_CHECK( loc_str != NULL, XLAL_EINVAL );
546 XLAL_CHECK( tmpl_str != NULL, XLAL_EINVAL );
547 XLAL_CHECK( param_tol_mism > 0, XLAL_EINVAL );
548 XLAL_CHECK( metric != NULL, XLAL_EFAULT );
549 XLAL_CHECK( rssky_transf != NULL, XLAL_EFAULT );
550 XLAL_CHECK( phys_1 != NULL, XLAL_EFAULT );
551 XLAL_CHECK( phys_2 != NULL, XLAL_EFAULT );
552
553 // Local copies of physical coordinates
554 const PulsarDopplerParams phys_orig[2] = { *phys_1, *phys_2 };
555 PulsarDopplerParams phys[2] = { phys_orig[0], phys_orig[1] };
556
557 // Round physical coordinates
558 for ( size_t i = 0; i < 2; ++i ) {
559 phys[i].Alpha = round_to_dp_sf( phys[i].Alpha, round_param_to_dp, round_param_to_sf );
560 phys[i].Delta = round_to_dp_sf( phys[i].Delta, round_param_to_dp, round_param_to_sf );
561 phys[i].fkdot[0] = round_to_dp_sf( phys[i].fkdot[0], round_param_to_dp, round_param_to_sf );
562 phys[i].fkdot[1] = round_to_dp_sf( phys[i].fkdot[1], round_param_to_dp, round_param_to_sf );
563 }
564
565 // Transform physical point to reduced supersky coordinates
566 double rssky_array[2][metric->size1];
567 gsl_vector_view rssky_view[2] = { gsl_vector_view_array( rssky_array[0], metric->size1 ), gsl_vector_view_array( rssky_array[1], metric->size1 ) };
568 gsl_vector *const rssky[2] = { &rssky_view[0].vector, &rssky_view[1].vector };
569 for ( size_t i = 0; i < 2; ++i ) {
570 XLAL_CHECK( XLALConvertPhysicalToSuperskyPoint( rssky[i], &phys[i], rssky_transf ) == XLAL_SUCCESS, XLAL_EFUNC );
571 }
572
573
574 // Store difference between reduced supersky coordinates in 'u'
575 double u_array[metric->size1];
576 gsl_vector_view u_view = gsl_vector_view_array( u_array, metric->size1 );
577 gsl_vector *const u = &u_view.vector;
578 gsl_vector_memcpy( u, rssky[0] );
579 gsl_vector_sub( u, rssky[1] );
580
581 // Multiply 'u' by metric, storing result in 'v'
582 double v_array[metric->size1];
583 gsl_vector_view v_view = gsl_vector_view_array( v_array, metric->size1 );
584 gsl_vector *const v = &v_view.vector;
585 gsl_blas_dsymv( CblasUpper, 1.0, metric, u, 0.0, v );
586
587 // Compute mismatch and compare to tolerance
588 REAL8 mism = 0;
589 gsl_blas_ddot( u, v, &mism );
590
591 // If mismatch is above tolerance, print error message
592 if ( mism > param_tol_mism ) {
593 *equal = 0;
594 XLALPrintInfo( "%s: at %s, mismatch between %s template parameters exceeds tolerance: %g > %g\n", __func__, loc_str, tmpl_str, mism, param_tol_mism );
595 XLALPrintInfo( "%s: serial 1 = %"LAL_UINT8_FORMAT"\n", __func__, serial_1 );
596 XLALPrintInfo( "%s: serial 2 = %"LAL_UINT8_FORMAT"\n", __func__, serial_2 );
597 for ( size_t i = 0; i < 2; ++i ) {
598 XLALPrintInfo( "%s: physical %zu = {%.15g,%.15g,%.15g,%.15g}\n", __func__, i + 1,
599 phys_orig[i].Alpha, phys_orig[i].Delta, phys_orig[i].fkdot[0], phys_orig[i].fkdot[1] );
600 if ( round_param_to_dp > 0 ) {
601 XLALPrintInfo( "%s: physical %zu = {%.15g,%.15g,%.15g,%.15g} rounded to %u d.p., %u s.f.\n", __func__, i + 1,
602 phys[i].Alpha, phys[i].Delta, phys[i].fkdot[0], phys[i].fkdot[1],
603 round_param_to_dp, round_param_to_sf );
604 }
605 }
606 for ( size_t i = 0; i < 2; ++i ) {
607 XLALPrintInfo( "%s: reduced supersky %zu = ", __func__, i + 1 );
608 for ( size_t j = 0; j < rssky[i]->size; ++j ) {
609 XLALPrintInfo( "%c%.15g", j == 0 ? '{' : ',', gsl_vector_get( rssky[i], j ) );
610 }
611 XLALPrintInfo( "}\n" );
612 }
613 XLALPrintInfo( "%s: reduced supersky diff = ", __func__ );
614 for ( size_t j = 0; j < u->size; ++j ) {
615 XLALPrintInfo( "%c%.15g", j == 0 ? '{' : ',', gsl_vector_get( u, j ) );
616 }
617 XLALPrintInfo( "}\n" );
618 XLALPrintInfo( "%s: metric dot = ", __func__ );
619 for ( size_t j = 0; j < u->size; ++j ) {
620 XLALPrintInfo( "%c%.15g", j == 0 ? '{' : ',', gsl_vector_get( u, j ) * gsl_vector_get( v, j ) );
621 }
622 XLALPrintInfo( "}\n" );
623 }
624
625 return XLAL_SUCCESS;
626
627}
628
629///
630/// Compare two vectors of results
631///
633 BOOLEAN *equal,
634 const VectorComparison *result_tol,
635 const REAL4Vector *res_1,
636 const REAL4Vector *res_2,
637 const UINT4 r1,
638 const UINT4 r2
639)
640{
641 XLAL_CHECK( r1 <= res_1->length, XLAL_EINVAL );
642 XLAL_CHECK( r2 <= res_2->length, XLAL_EINVAL );
643 const REAL4Vector res_1_n = { .length = r1, .data = res_1->data };
644 const REAL4Vector res_2_n = { .length = r2, .data = res_2->data };
645 VectorComparison XLAL_INIT_DECL( result_diff );
646 int errnum = 0;
647 XLAL_TRY( XLALCompareREAL4Vectors( &result_diff, &res_1_n, &res_2_n, result_tol ), errnum );
648 if ( errnum & XLAL_ETOL ) {
649 *equal = 0;
650 } else if ( errnum != 0 ) {
652 }
653 return XLAL_SUCCESS;
654}
655
656///
657/// Create results toplist
658///
659WeaveResultsToplist *XLALWeaveResultsToplistCreate(
660 const size_t nspins,
661 WeaveStatisticsParams *statistics_params,
662 const char *stat_name,
663 const char *stat_desc,
664 const UINT4 toplist_limit,
665 WeaveResultsToplistRankingStats toplist_rank_stats_fcn,
666 WeaveResultsToplistItemGetRankStat toplist_item_get_rank_stat_fcn,
667 WeaveResultsToplistItemSetRankStat toplist_item_set_rank_stat_fcn
668)
669{
670
671 // Check input
672 XLAL_CHECK_NULL( stat_name != NULL, XLAL_EFAULT );
673 XLAL_CHECK_NULL( stat_desc != NULL, XLAL_EFAULT );
674 XLAL_CHECK_NULL( statistics_params != NULL, XLAL_EFAULT );
675
676 // Allocate memory
677 WeaveResultsToplist *toplist = XLALCalloc( 1, sizeof( *toplist ) );
678 XLAL_CHECK_NULL( toplist != NULL, XLAL_ENOMEM );
679
680 // Set fields
681 toplist->nspins = nspins;
682 toplist->stat_name = stat_name;
683 toplist->stat_desc = stat_desc;
684 toplist->rank_stats_fcn = toplist_rank_stats_fcn;
685 toplist->item_get_rank_stat_fcn = toplist_item_get_rank_stat_fcn;
686 toplist->item_set_rank_stat_fcn = toplist_item_set_rank_stat_fcn;
687 toplist->statistics_params = statistics_params;
688
689 // Create heap which ranks toplist items
690 toplist->heap = XLALHeapCreate2( ( LALHeapDtorFcn ) toplist_item_destroy, toplist_limit, +1, toplist_item_compare, toplist_item_get_rank_stat_fcn );
691 XLAL_CHECK_NULL( toplist->heap != NULL, XLAL_EFUNC );
692
693 return toplist;
694
695}
696
697///
698/// Free results toplist
699///
701 WeaveResultsToplist *toplist
702)
703{
704 if ( toplist != NULL ) {
705 XLALDestroyUINT4Vector( toplist->maybe_add_freq_idxs );
706 XLALHeapDestroy( toplist->heap );
707 toplist_item_destroy( toplist->saved_item );
708 XLALFree( toplist );
709 }
710}
711
712///
713/// Add semicoherent results to toplist
714///
716 WeaveResultsToplist *toplist,
717 const WeaveSemiResults *semi_res,
718 const UINT4 semi_nfreqs
719)
720{
721 // Check input
722 XLAL_CHECK( toplist != NULL, XLAL_EFAULT );
723 XLAL_CHECK( semi_res != NULL, XLAL_EFAULT );
724
725 // Reallocate vector of indexes of toplist results which should be considered for addition
726 if ( toplist->maybe_add_freq_idxs == NULL || toplist->maybe_add_freq_idxs->length < semi_nfreqs ) {
727 toplist->maybe_add_freq_idxs = XLALResizeUINT4Vector( toplist->maybe_add_freq_idxs, semi_nfreqs );
728 XLAL_CHECK( toplist->maybe_add_freq_idxs != NULL, XLAL_EFUNC );
729 }
730
731 // Get pointer to array of ranking statistics
732 const REAL4 *toplist_rank_stats = toplist->rank_stats_fcn( semi_res );
733
734 // Get ranking statistic of heap root (or -infinity if heap is not yet full)
735 const int heap_full = XLALHeapIsFull( toplist->heap );
736 XLAL_CHECK( heap_full >= 0, XLAL_EFUNC );
737 const REAL4 heap_root_rank_stat = heap_full ? toplist->item_get_rank_stat_fcn( XLALHeapRoot( toplist->heap ) ) : GSL_NEGINF;
738
739 // Find the indexes of the semicoherent results whose ranking statistic equals or exceeds
740 // that of the heap root; only select these results for possible insertion into the toplist
741 UINT4 n_maybe_add = 0;
742 XLAL_CHECK( XLALVectorFindScalarLessEqualREAL4( &n_maybe_add, toplist->maybe_add_freq_idxs->data, heap_root_rank_stat, toplist_rank_stats, semi_nfreqs ) == XLAL_SUCCESS, XLAL_EFUNC );
743
744 // Whether we output per-segment template coordinates is currently tied to output of any per-segment statistics
745 const WeaveStatisticsParams *params = toplist->statistics_params;
746 WeaveStatisticType per_seg_coords = params->statistics_to_output[0] & ( WEAVE_STATISTIC_COH2F | WEAVE_STATISTIC_COH2F_DET );
747
748 // Iterate over semicoherent results which have been selected for possible toplist insertion
749 for ( UINT4 idx = 0; idx < n_maybe_add; ++idx ) {
750 const UINT4 freq_idx = toplist->maybe_add_freq_idxs->data[idx];
751
752 // Create a new toplist item if needed
753 if ( toplist->saved_item == NULL ) {
754 toplist->saved_item = toplist_item_create( toplist );
755 XLAL_CHECK( toplist->saved_item != NULL, XLAL_ENOMEM );
756 }
757 WeaveResultsToplistItem *item = toplist->saved_item;
758
759 // Set ranking statistic of toplist item
760 toplist->item_set_rank_stat_fcn( item, toplist_rank_stats[freq_idx] );
761
762 // Possibly add toplist item to heap
763 XLAL_CHECK( XLALHeapAdd( toplist->heap, ( void ** ) &toplist->saved_item ) == XLAL_SUCCESS, XLAL_EFUNC );
764
765 // Skip remainder of loop if toplist item was not added to heap
766 if ( item == toplist->saved_item ) {
767 continue;
768 }
769
770 // Set serial number of template
771 item->serial = ++toplist->serial;
772
773 // Set all semicoherent template parameters
774 item->semi_alpha = semi_res->semi_phys.Alpha;
775 item->semi_delta = semi_res->semi_phys.Delta;
776 item->semi_fkdot[0] = semi_res->semi_phys.fkdot[0] + freq_idx * semi_res->dfreq;
777 for ( size_t k = 1; k <= toplist->nspins; ++k ) {
778 item->semi_fkdot[k] = semi_res->semi_phys.fkdot[k];
779 }
780
781 // Set all coherent template parameters if outputting per-segment statistics
782 if ( per_seg_coords ) {
783 for ( size_t j = 0; j < semi_res->nsegments; ++j ) {
784 item->coh_alpha[j] = semi_res->coh_phys[j].Alpha;
785 item->coh_delta[j] = semi_res->coh_phys[j].Delta;
786 item->coh_fkdot[0][j] = semi_res->coh_phys[j].fkdot[0] + freq_idx * semi_res->dfreq;
787 for ( size_t k = 1; k <= toplist->nspins; ++k ) {
788 item->coh_fkdot[k][j] = semi_res->coh_phys[j].fkdot[k];
789 }
790 }
791 }
792
793 // Skip remainder of loop if simulating search
794 if ( semi_res->simulation_level & WEAVE_SIMULATE ) {
795 continue;
796 }
797
798 //
799 // Copy all 'mainloop_statistics_to_keep' statistic values, as they will be needed either 1) for output or 2) computing remaining completion-loop statistics
800 //
801 WeaveStatisticType stats_to_keep = params->mainloop_statistics_to_keep;
802
803 if ( stats_to_keep & WEAVE_STATISTIC_COH2F ) {
804 XLAL_CHECK( XLALWeaveSemiCoh2FExtract( item->stage[0].coh2F, semi_res, freq_idx ) == XLAL_SUCCESS, XLAL_EFUNC );
805 }
806 if ( stats_to_keep & WEAVE_STATISTIC_COH2F_DET ) {
807 for ( size_t i = 0; i < semi_res->ndetectors; ++i ) {
808 for ( size_t j = 0; j < semi_res->nsegments; ++j ) {
809 if ( semi_res->coh2F_det[i][j] != NULL ) {
810 item->stage[0].coh2F_det[i][j] = semi_res->coh2F_det[i][j][freq_idx];
811 } else {
812 // There is not per-detector F-statistic for this segment, usually because this segment contains
813 // no data from this detector. In this case we output a clearly invalid F-statistic value.
814 item->stage[0].coh2F_det[i][j] = NAN;
815 }
816 }
817 }
818 }
819
820 if ( stats_to_keep & WEAVE_STATISTIC_MAX2F ) {
821 item->stage[0].max2F = semi_res->max2F->data[freq_idx];
822 }
823 if ( stats_to_keep & WEAVE_STATISTIC_MAX2F_DET ) {
824 for ( size_t i = 0; i < semi_res->ndetectors; ++i ) {
825 item->stage[0].max2F_det[i] = semi_res->max2F_det[i]->data[freq_idx];
826 }
827 }
828
829 if ( stats_to_keep & WEAVE_STATISTIC_SUM2F ) {
830 item->stage[0].sum2F = semi_res->sum2F->data[freq_idx];
831 }
832 if ( stats_to_keep & WEAVE_STATISTIC_SUM2F_DET ) {
833 for ( size_t i = 0; i < semi_res->ndetectors; ++i ) {
834 item->stage[0].sum2F_det[i] = semi_res->sum2F_det[i]->data[freq_idx];
835 }
836 }
837
838 if ( stats_to_keep & WEAVE_STATISTIC_MEAN2F ) {
839 item->stage[0].mean2F = semi_res->mean2F->data[freq_idx];
840 }
841
842 if ( stats_to_keep & WEAVE_STATISTIC_BSGL ) {
843 item->stage[0].log10BSGL = semi_res->log10BSGL->data[freq_idx];
844 }
845
846 if ( stats_to_keep & WEAVE_STATISTIC_BSGLtL ) {
847 item->stage[0].log10BSGLtL = semi_res->log10BSGLtL->data[freq_idx];
848 }
849
850 if ( stats_to_keep & WEAVE_STATISTIC_BtSGLtL ) {
851 item->stage[0].log10BtSGLtL = semi_res->log10BtSGLtL->data[freq_idx];
852 }
853
854 }
855
856 return XLAL_SUCCESS;
857
858}
859
860///
861/// Compute all missing 'extra' (non-toplist-ranking) statistics for all toplist entries
862///
864 WeaveResultsToplist *toplist
865)
866{
867 // Check input
868 XLAL_CHECK( toplist != NULL, XLAL_EFAULT );
869
870 // Compute all completion-loop statistics on toplist items
871 XLAL_CHECK( XLALHeapModify( toplist->heap, toplist_fill_completionloop_stats, toplist->statistics_params ) == XLAL_SUCCESS, XLAL_EFUNC );
872
873 return XLAL_SUCCESS;
874
875}
876
877///
878/// Write results toplist to a FITS file
879///
881 FITSFile *file,
882 const WeaveResultsToplist *toplist
883)
884{
885
886 // Check input
887 XLAL_CHECK( file != NULL, XLAL_EFAULT );
888 XLAL_CHECK( toplist != NULL, XLAL_EFAULT );
889
890 // Format name and description of statistic
891 char name[256];
892 snprintf( name, sizeof( name ), "%s_toplist", toplist->stat_name );
893 char desc[256];
894 snprintf( desc, sizeof( desc ), "toplist ranked by %s", toplist->stat_desc );
895
896 // Open FITS table for writing and initialise
899
900 // Write all heap items to FITS table
902
903 // Write current toplist item serial
904 XLAL_CHECK( XLALFITSHeaderWriteUINT8( file, "serial", toplist->serial, "item serial" ) == XLAL_SUCCESS, XLAL_EFUNC );
905
906 return XLAL_SUCCESS;
907
908}
909
910///
911/// Read results from a FITS file and append to existing results toplist
912///
914 FITSFile *file,
915 WeaveResultsToplist *toplist
916)
917{
918
919 // Check input
920 XLAL_CHECK( file != NULL, XLAL_EFAULT );
921 XLAL_CHECK( toplist != NULL, XLAL_EFAULT );
922
923 // Format name of statistic
924 char name[256];
925 snprintf( name, sizeof( name ), "%s_toplist", toplist->stat_name );
926
927 // Open FITS table for reading and initialise
928 UINT8 nrows = 0;
931
932 // Read all items from FITS table
933 while ( nrows > 0 ) {
934
935 // Create a new toplist item if needed
936 if ( toplist->saved_item == NULL ) {
937 toplist->saved_item = toplist_item_create( toplist );
938 XLAL_CHECK( toplist->saved_item != NULL, XLAL_ENOMEM );
939 }
940
941 // Read item from FITS table
942 XLAL_CHECK( XLALFITSTableReadRow( file, toplist->saved_item, &nrows ) == XLAL_SUCCESS, XLAL_EFUNC );
943
944 // Save highest serial in toplist
945 if ( toplist->serial < toplist->saved_item->serial ) {
946 toplist->serial = toplist->saved_item->serial;
947 }
948
949 // Add item to heap
950 XLAL_CHECK( XLALHeapAdd( toplist->heap, ( void ** ) &toplist->saved_item ) == XLAL_SUCCESS, XLAL_EFUNC );
951
952 }
953
954 // Write current toplist item serial
955 XLAL_CHECK( XLALFITSHeaderReadUINT8( file, "serial", &toplist->serial ) == XLAL_SUCCESS, XLAL_EFUNC );
956
957 return XLAL_SUCCESS;
958
959}
960
961///
962/// Compare two results toplists and return whether they are equal
963///
965 BOOLEAN *equal,
966 const WeaveSetupData *setup,
967 const BOOLEAN sort_by_semi_phys,
968 const UINT4 round_param_to_dp,
969 const UINT4 round_param_to_sf,
970 const REAL8 unmatched_item_tol,
971 const REAL8 param_tol_mism,
972 const VectorComparison *result_tol,
973 const UINT4 toplist_compare_limit,
974 const WeaveResultsToplist *toplist_1,
975 const WeaveResultsToplist *toplist_2
976)
977{
978
979 // Check input
980 XLAL_CHECK( equal != NULL, XLAL_EFAULT );
981 XLAL_CHECK( setup != NULL, XLAL_EFAULT );
982 XLAL_CHECK( unmatched_item_tol >= 0, XLAL_EINVAL );
983 XLAL_CHECK( param_tol_mism >= 0, XLAL_EINVAL );
984 XLAL_CHECK( result_tol != NULL, XLAL_EFAULT );
985 XLAL_CHECK( toplist_1 != NULL, XLAL_EFAULT );
986 XLAL_CHECK( toplist_2 != NULL, XLAL_EFAULT );
987 XLAL_CHECK( strcmp( toplist_1->stat_name, toplist_2->stat_name ) == 0, XLAL_EINVAL );
988 XLAL_CHECK( strcmp( toplist_1->stat_desc, toplist_2->stat_desc ) == 0, XLAL_EINVAL );
989
990 const WeaveResultsToplist *toplist = toplist_1;
991 const WeaveStatisticsParams *params = toplist->statistics_params;
992
993 // Results toplists are assumed equal until we find otherwise
994 *equal = 1;
995
996 // Compare lengths of heaps
997 const size_t n = XLALHeapSize( toplist_1->heap );
998 {
999 const size_t n_2 = XLALHeapSize( toplist_2->heap );
1000 if ( n != n_2 ) {
1001 *equal = 0;
1002 XLALPrintInfo( "%s: unequal size %s toplists: %zu != %zu\n", __func__, toplist->stat_desc, n, n_2 );
1003 return XLAL_SUCCESS;
1004 }
1005 }
1006
1007 // Get lists of toplist items
1008 const WeaveResultsToplistItem **items_1 = ( const WeaveResultsToplistItem ** ) XLALHeapElements( toplist_1->heap );
1009 XLAL_CHECK( items_1 != NULL, XLAL_EFUNC );
1010 const WeaveResultsToplistItem **items_2 = ( const WeaveResultsToplistItem ** ) XLALHeapElements( toplist_2->heap );
1011 XLAL_CHECK( items_2 != NULL, XLAL_EFUNC );
1012
1013 // Build list of matched items
1014 const WeaveResultsToplistItem **matched_1 = XLALCalloc( n, sizeof( *matched_1 ) );
1015 XLAL_CHECK( matched_1 != NULL, XLAL_ENOMEM );
1016 const WeaveResultsToplistItem **matched_2 = XLALCalloc( n, sizeof( *matched_2 ) );
1017 XLAL_CHECK( matched_2 != NULL, XLAL_ENOMEM );
1018 size_t matched_n = 0;
1019 if ( sort_by_semi_phys ) {
1020
1021 // Sort toplist items by physical coordinates of semicoherent template
1022 qsort( items_1, n, sizeof( *items_1 ), toplist_item_sort_by_semi_phys );
1023 qsort( items_2, n, sizeof( *items_2 ), toplist_item_sort_by_semi_phys );
1024
1025 // Assume this order is matched
1026 memcpy( matched_1, items_1, n * sizeof( *matched_1 ) );
1027 memcpy( matched_2, items_2, n * sizeof( *matched_2 ) );
1028 matched_n = n;
1029
1030 } else {
1031
1032 // Sort toplist items by serial number of template
1033 qsort( items_1, n, sizeof( *items_1 ), toplist_item_sort_by_serial );
1034 qsort( items_2, n, sizeof( *items_2 ), toplist_item_sort_by_serial );
1035
1036 // Match up items by serial number
1037 size_t n_1 = 0, n_2 = 0;
1038 while ( n_1 < n && n_2 < n ) {
1039 if ( items_1[n_1]->serial < items_2[n_2]->serial ) {
1040 ++n_1;
1041 } else if ( items_2[n_2]->serial < items_1[n_1]->serial ) {
1042 ++n_2;
1043 } else {
1044 matched_1[matched_n] = items_1[n_1++];
1045 matched_2[matched_n] = items_2[n_2++];
1046 ++matched_n;
1047 }
1048 }
1049
1050 }
1051
1052 // Limit number of items to match
1053 if ( toplist_compare_limit > 0 && matched_n > toplist_compare_limit ) {
1054 matched_n = toplist_compare_limit;
1055 } else {
1056
1057 // Check that number of unmatched toplist items does not exceed tolerance
1058 const long int unmatched_n = n - matched_n;
1059 const long int unmatched_item_tol_n = lrint( unmatched_item_tol * n );
1060 if ( unmatched_n > unmatched_item_tol_n ) {
1061 *equal = 0;
1062 XLALPrintInfo( "%s: number of unmatched %s toplist items %lu exceeds tolerance %0.3e * %zu = %lu\n", __func__, toplist->stat_desc, unmatched_n, unmatched_item_tol, n, unmatched_item_tol_n );
1063 }
1064
1065 }
1066
1067 // Allocate vectors for storing results for comparison with compare_vectors()
1068 REAL4Vector *res_1 = XLALCreateREAL4Vector( matched_n * params->detectors->length * params->nsegments );
1069 XLAL_CHECK( res_1 != NULL, XLAL_EFUNC );
1070 REAL4Vector *res_2 = XLALCreateREAL4Vector( matched_n * params->detectors->length * params->nsegments );
1071 XLAL_CHECK( res_2 != NULL, XLAL_EFUNC );
1072
1073 // Print toplists being compared and number of matched items
1074 XLALPrintInfo( "%s: comparing %zu matched items from toplists ranked by %s ...\n", __func__, matched_n, toplist->stat_desc );
1075
1076 while ( *equal ) { // So we can use 'break' to skip comparisons on failure
1077
1078 // Compare semicoherent template parameters
1079 if ( ( *equal ) && ( param_tol_mism > 0 ) ) {
1080 for ( size_t i = 0; i < matched_n; ++i ) {
1081 char loc_str[256];
1082 snprintf( loc_str, sizeof( loc_str ), "toplist item %zu", i );
1083 const UINT8 serial_1 = matched_1[i]->serial;
1084 const UINT8 serial_2 = matched_2[i]->serial;
1085 PulsarDopplerParams XLAL_INIT_DECL( semi_phys_1 );
1086 PulsarDopplerParams XLAL_INIT_DECL( semi_phys_2 );
1087 semi_phys_1.Alpha = matched_1[i]->semi_alpha;
1088 semi_phys_2.Alpha = matched_2[i]->semi_alpha;
1089 semi_phys_1.Delta = matched_1[i]->semi_delta;
1090 semi_phys_2.Delta = matched_2[i]->semi_delta;
1091 for ( size_t k = 0; k <= toplist->nspins; ++k ) {
1092 semi_phys_1.fkdot[k] = matched_1[i]->semi_fkdot[k];
1093 semi_phys_2.fkdot[k] = matched_2[i]->semi_fkdot[k];
1094 };
1096 loc_str, "semicoherent",
1097 round_param_to_dp, round_param_to_sf, param_tol_mism,
1098 setup->metrics->semi_rssky_metric, setup->metrics->semi_rssky_transf,
1099 serial_1, serial_2,
1100 &semi_phys_1, &semi_phys_2
1101 ) == XLAL_SUCCESS, XLAL_EFUNC );
1102 if ( !*equal ) {
1103 break;
1104 }
1105 }
1106 if ( !*equal ) {
1107 break;
1108 }
1109 }
1110
1111 // Compare coherent template parameters
1112 if ( ( *equal ) && ( param_tol_mism > 0 ) && ( params->statistics_to_output[0] & ( WEAVE_STATISTIC_COH2F | WEAVE_STATISTIC_COH2F_DET ) ) ) {
1113 for ( size_t i = 0; i < matched_n; ++i ) {
1114 for ( size_t j = 0; j < params->nsegments; ++j ) {
1115 char loc_str[256];
1116 snprintf( loc_str, sizeof( loc_str ), "toplist item %zu, segment %zu", i, j );
1117 const UINT8 serial_1 = matched_1[i]->serial;
1118 const UINT8 serial_2 = matched_2[i]->serial;
1119 PulsarDopplerParams XLAL_INIT_DECL( coh_phys_1 );
1120 PulsarDopplerParams XLAL_INIT_DECL( coh_phys_2 );
1121 coh_phys_1.Alpha = matched_1[i]->coh_alpha[j];
1122 coh_phys_2.Alpha = matched_2[i]->coh_alpha[j];
1123 coh_phys_1.Delta = matched_1[i]->coh_delta[j];
1124 coh_phys_2.Delta = matched_2[i]->coh_delta[j];
1125 for ( size_t k = 0; k <= toplist->nspins; ++k ) {
1126 coh_phys_1.fkdot[k] = matched_1[i]->coh_fkdot[k][j];
1127 coh_phys_2.fkdot[k] = matched_2[i]->coh_fkdot[k][j];
1128 };
1130 loc_str, "coherent",
1131 round_param_to_dp, round_param_to_sf, param_tol_mism,
1132 setup->metrics->coh_rssky_metric[j], setup->metrics->coh_rssky_transf[j],
1133 serial_1, serial_2,
1134 &coh_phys_1, &coh_phys_2
1135 ) == XLAL_SUCCESS, XLAL_EFUNC );
1136 }
1137 if ( !*equal ) {
1138 break;
1139 }
1140 }
1141 if ( !*equal ) {
1142 break;
1143 }
1144 }
1145
1146 // Compare statistics values from both stages ('stage 0' = main search, 'stage 1' = recalculation stage ('recalc'))
1147 for ( UINT4 istage = 0; istage < 2; ++ istage ) {
1148 WeaveStatisticType stats_to_output = params->statistics_to_output[istage];
1149
1150 // Compare mean multi-detector F-statistics
1151 if ( stats_to_output & WEAVE_STATISTIC_MEAN2F ) {
1152 XLALPrintInfo( "%s: comparing mean multi-detector F-statistics ...\n", __func__ );
1153 UINT4 r1 = 0, r2 = 0;
1154 for ( size_t i = 0; i < matched_n; ++i ) {
1155 res_1->data[r1++] = matched_1[i]->stage[istage].mean2F;
1156 res_2->data[r2++] = matched_2[i]->stage[istage].mean2F;
1157 }
1158 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1159 }
1160
1161 // Compare mean per-detector F-statistic
1162 if ( stats_to_output & WEAVE_STATISTIC_MEAN2F_DET ) {
1163 XLALPrintInfo( "%s: comparing mean per-detector F-statistics ...\n", __func__ );
1164 UINT4 r1 = 0, r2 = 0;
1165 for ( size_t k = 0; k < params->detectors->length; ++k ) {
1166 for ( size_t i = 0; i < matched_n; ++i ) {
1167 res_1->data[r1++] = matched_1[i]->stage[istage].mean2F_det[k];
1168 res_2->data[r2++] = matched_2[i]->stage[istage].mean2F_det[k];
1169 }
1170 }
1171 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1172 }
1173
1174 // Compare per-segment coherent multi-detector F-statistics
1175 if ( stats_to_output & WEAVE_STATISTIC_COH2F ) {
1176 XLALPrintInfo( "%s: comparing coherent multi-detector F-statistics ...\n", __func__ );
1177 UINT4 r1 = 0, r2 = 0;
1178 for ( size_t j = 0; j < params->nsegments; ++j ) {
1179 for ( size_t i = 0; i < matched_n; ++i ) {
1180 res_1->data[r1++] = matched_1[i]->stage[istage].coh2F[j];
1181 res_2->data[r2++] = matched_2[i]->stage[istage].coh2F[j];
1182 }
1183 }
1184 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1185 }
1186
1187 // Compare per-segment per-detector F-statistics
1188 if ( stats_to_output & WEAVE_STATISTIC_COH2F_DET ) {
1189 XLALPrintInfo( "%s: comparing per-segment per-detector F-statistics ...\n", __func__ );
1190 UINT4 r1 = 0, r2 = 0;
1191 for ( size_t j = 0; j < params->nsegments; ++j ) {
1192 for ( size_t k = 0; k < params->detectors->length; ++k ) {
1193 if ( isfinite( items_1[0]->stage[istage].coh2F_det[k][j] ) || isfinite( items_2[0]->stage[istage].coh2F_det[k][j] ) ) {
1194 for ( size_t i = 0; i < matched_n; ++i ) {
1195 res_1->data[r1++] = matched_1[i]->stage[istage].coh2F_det[k][j];
1196 res_2->data[r2++] = matched_2[i]->stage[istage].coh2F_det[k][j];
1197 }
1198 }
1199 }
1200 }
1201 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1202 }
1203
1204 // Compare segment-max multi-detector F-statistics
1205 if ( stats_to_output & WEAVE_STATISTIC_MAX2F ) {
1206 XLALPrintInfo( "%s: comparing max multi-detector F-statistics ...\n", __func__ );
1207 UINT4 r1 = 0, r2 = 0;
1208 for ( size_t i = 0; i < matched_n; ++i ) {
1209 res_1->data[r1++] = matched_1[i]->stage[istage].max2F;
1210 res_2->data[r2++] = matched_2[i]->stage[istage].max2F;
1211 }
1212 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1213 }
1214
1215 // Compare segment-max per-detector F-statistic
1216 if ( stats_to_output & WEAVE_STATISTIC_MAX2F_DET ) {
1217 XLALPrintInfo( "%s: comparing max per-detector F-statistics ...\n", __func__ );
1218 UINT4 r1 = 0, r2 = 0;
1219 for ( size_t k = 0; k < params->detectors->length; ++k ) {
1220 for ( size_t i = 0; i < matched_n; ++i ) {
1221 res_1->data[r1++] = matched_1[i]->stage[istage].max2F_det[k];
1222 res_2->data[r2++] = matched_2[i]->stage[istage].max2F_det[k];
1223 }
1224 }
1225 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1226 }
1227
1228 // Compare summed multi-detector F-statistics
1229 if ( stats_to_output & WEAVE_STATISTIC_SUM2F ) {
1230 XLALPrintInfo( "%s: comparing sum multi-detector F-statistics ...\n", __func__ );
1231 UINT4 r1 = 0, r2 = 0;
1232 for ( size_t i = 0; i < matched_n; ++i ) {
1233 res_1->data[r1++] = matched_1[i]->stage[istage].sum2F;
1234 res_2->data[r2++] = matched_2[i]->stage[istage].sum2F;
1235 }
1236 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1237 }
1238
1239 // Compare sum per-detector F-statistic
1240 if ( stats_to_output & WEAVE_STATISTIC_SUM2F_DET ) {
1241 XLALPrintInfo( "%s: comparing sum per-detector F-statistics ...\n", __func__ );
1242 UINT4 r1 = 0, r2 = 0;
1243 for ( size_t k = 0; k < params->detectors->length; ++k ) {
1244 for ( size_t i = 0; i < matched_n; ++i ) {
1245 res_1->data[r1++] = matched_1[i]->stage[istage].sum2F_det[k];
1246 res_2->data[r2++] = matched_2[i]->stage[istage].sum2F_det[k];
1247 }
1248 }
1249 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1250 }
1251
1252 // Compare line-robust BSGL statistic
1253 if ( stats_to_output & WEAVE_STATISTIC_BSGL ) {
1254 XLALPrintInfo( "%s: comparing line-robust B_S/GL statistic ...\n", __func__ );
1255 UINT4 r1 = 0, r2 = 0;
1256 for ( size_t i = 0; i < matched_n; ++i ) {
1257 res_1->data[r1++] = matched_1[i]->stage[istage].log10BSGL;
1258 res_2->data[r2++] = matched_2[i]->stage[istage].log10BSGL;
1259 }
1260 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1261 }
1262
1263 // Compare transient line-robust BSGLtL statistic
1264 if ( stats_to_output & WEAVE_STATISTIC_BSGLtL ) {
1265 XLALPrintInfo( "%s: comparing transient line-robust B_S/GLtL statistic ...\n", __func__ );
1266 UINT4 r1 = 0, r2 = 0;
1267 for ( size_t i = 0; i < matched_n; ++i ) {
1268 res_1->data[r1++] = matched_1[i]->stage[istage].log10BSGLtL;
1269 res_2->data[r2++] = matched_2[i]->stage[istage].log10BSGLtL;
1270 }
1271 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1272 }
1273
1274 // Compare transient signal line-robust BtSGLtL statistic
1275 if ( stats_to_output & WEAVE_STATISTIC_BtSGLtL ) {
1276 XLALPrintInfo( "%s: comparing transient signal line-robust B_tS/GLtL statistic ...\n", __func__ );
1277 UINT4 r1 = 0, r2 = 0;
1278 for ( size_t i = 0; i < matched_n; ++i ) {
1279 res_1->data[r1++] = matched_1[i]->stage[istage].log10BtSGLtL;
1280 res_2->data[r2++] = matched_2[i]->stage[istage].log10BtSGLtL;
1281 }
1282 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1283 }
1284
1285 // Compare 'Hough' multi-detector line statistics
1286 if ( stats_to_output & WEAVE_STATISTIC_NCOUNT ) {
1287 XLALPrintInfo( "%s: comparing 'Hough' multi-detector number count statistic ...\n", __func__ );
1288 UINT4 r1 = 0, r2 = 0;
1289 for ( size_t i = 0; i < matched_n; ++i ) {
1290 res_1->data[r1++] = matched_1[i]->stage[istage].ncount;
1291 res_2->data[r2++] = matched_2[i]->stage[istage].ncount;
1292 }
1293 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1294 }
1295
1296 // Compare 'Hough' per-detector line statistics
1297 if ( stats_to_output & WEAVE_STATISTIC_NCOUNT_DET ) {
1298 XLALPrintInfo( "%s: comparing 'Hough' per-detector number-count statistic ...\n", __func__ );
1299 UINT4 r1 = 0, r2 = 0;
1300 for ( size_t k = 0; k < params->detectors->length; ++k ) {
1301 for ( size_t i = 0; i < matched_n; ++i ) {
1302 res_1->data[r1++] = matched_1[i]->stage[istage].ncount_det[k];
1303 res_2->data[r2++] = matched_2[i]->stage[istage].ncount_det[k];
1304 }
1305 }
1306 XLAL_CHECK( compare_vectors( equal, result_tol, res_1, res_2, r1, r2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1307 }
1308
1309 if ( !*equal ) {
1310 break;
1311 }
1312
1313 } // istage
1314
1315 // All comparisons are done
1316 break;
1317
1318 } // while ( *equal )
1319
1320 // Cleanup
1321 XLALFree( items_1 );
1322 XLALFree( items_2 );
1323 XLALFree( matched_1 );
1324 XLALFree( matched_2 );
1325 XLALDestroyREAL4Vector( res_1 );
1326 XLALDestroyREAL4Vector( res_2 );
1327
1328 return XLAL_SUCCESS;
1329
1330}
1331
1332// Local Variables:
1333// c-file-style: "linux"
1334// c-basic-offset: 2
1335// End:
int XLALWeaveCohResultsExtract(REAL4Vector **coh2F, REAL4Vector *coh2F_det[PULSAR_MAX_DETECTORS], BOOLEAN *have_coh2F_det, WeaveCohResults *coh_res, const WeaveCohInput *coh_input)
Simple API function to extract pointers to 2F results from WeaveCohResults.
int XLALWeaveCohResultsCompute(WeaveCohResults **coh_res, WeaveCohInput *coh_input, const PulsarDopplerParams *coh_phys, const UINT4 coh_nfreqs, WeaveSearchTiming *tim)
Create and compute coherent results.
int XLALWeaveSemiCoh2FExtract(REAL4 *coh2F, const WeaveSemiResults *semi_res, const UINT4 freq_idx)
Extract 2F results from WeaveSemiResults; handles results stores in CUDA device memory.
Module which computes coherent and semicoherent results.
#define __func__
log an I/O error, i.e.
int XLALFITSTableWriteRow(FITSFile UNUSED *file, const void UNUSED *record)
Definition: FITSFileIO.c:2550
int XLALFITSTableReadRow(FITSFile UNUSED *file, void UNUSED *record, UINT8 UNUSED *rem_nrows)
Definition: FITSFileIO.c:2621
int XLALFITSTableOpenWrite(FITSFile UNUSED *file, const CHAR UNUSED *name, const CHAR UNUSED *comment)
Definition: FITSFileIO.c:2162
int XLALFITSHeaderWriteUINT8(FITSFile UNUSED *file, const CHAR UNUSED *key, const UINT8 UNUSED value, const CHAR UNUSED *comment)
Definition: FITSFileIO.c:827
int XLALFITSHeaderReadUINT8(FITSFile UNUSED *file, const CHAR UNUSED *key, UINT8 UNUSED *value)
Definition: FITSFileIO.c:863
int XLALFITSTableOpenRead(FITSFile UNUSED *file, const CHAR UNUSED *name, UINT8 UNUSED *nrows)
Definition: FITSFileIO.c:2198
int j
int k
stage
int XLALWeaveResultsToplistAdd(WeaveResultsToplist *toplist, const WeaveSemiResults *semi_res, const UINT4 semi_nfreqs)
Add semicoherent results to toplist.
static int toplist_fits_table_write_visitor(void *param, const void *x)
Visitor function for writing a toplist to a FITS table.
void XLALWeaveResultsToplistDestroy(WeaveResultsToplist *toplist)
Free results toplist.
static double round_to_dp_sf(double x, const UINT4 dp, const UINT4 sf)
Round a floating-point number 'x' to 'dp' decimal places (if > 0), then 'sf' significant figures (if ...
WeaveResultsToplist * XLALWeaveResultsToplistCreate(const size_t nspins, WeaveStatisticsParams *statistics_params, const char *stat_name, const char *stat_desc, const UINT4 toplist_limit, WeaveResultsToplistRankingStats toplist_rank_stats_fcn, WeaveResultsToplistItemGetRankStat toplist_item_get_rank_stat_fcn, WeaveResultsToplistItemSetRankStat toplist_item_set_rank_stat_fcn)
Create results toplist.
static int toplist_fits_table_init(FITSFile *file, const WeaveResultsToplist *toplist)
Initialise a FITS table for writing/reading a toplist.
static int toplist_item_sort_by_serial(const void *x, const void *y)
Sort toplist items by serial number of template.
int XLALWeaveResultsToplistReadAppend(FITSFile *file, WeaveResultsToplist *toplist)
Read results from a FITS file and append to existing results toplist.
static int toplist_item_compare(void *param, const void *x, const void *y)
Compare toplist items.
#define COMPARE_BY(x, y)
static void toplist_item_destroy(WeaveResultsToplistItem *item)
Destroy a toplist item.
int XLALWeaveResultsToplistCompare(BOOLEAN *equal, const WeaveSetupData *setup, const BOOLEAN sort_by_semi_phys, const UINT4 round_param_to_dp, const UINT4 round_param_to_sf, const REAL8 unmatched_item_tol, const REAL8 param_tol_mism, const VectorComparison *result_tol, const UINT4 toplist_compare_limit, const WeaveResultsToplist *toplist_1, const WeaveResultsToplist *toplist_2)
Compare two results toplists and return whether they are equal.
static int compare_vectors(BOOLEAN *equal, const VectorComparison *result_tol, const REAL4Vector *res_1, const REAL4Vector *res_2, const UINT4 r1, const UINT4 r2)
Compare two vectors of results.
static int compare_templates(BOOLEAN *equal, const char *loc_str, const char *tmpl_str, const UINT4 round_param_to_dp, const UINT4 round_param_to_sf, const REAL8 param_tol_mism, const gsl_matrix *metric, const SuperskyTransformData *rssky_transf, const UINT8 serial_1, const UINT8 serial_2, const PulsarDopplerParams *phys_1, const PulsarDopplerParams *phys_2)
Compute two template parameters.
int XLALWeaveResultsToplistWrite(FITSFile *file, const WeaveResultsToplist *toplist)
Write results toplist to a FITS file.
static WeaveResultsToplistItem * toplist_item_create(const WeaveResultsToplist *toplist)
Create a toplist item.
static int toplist_fill_completionloop_stats(void *param, void *x)
Function to update given toplist item with missing 'completion loop' statistics.
static int toplist_item_sort_by_semi_phys(const void *x, const void *y)
Sort toplist items by physical coordinates of semicoherent template.
int XLALWeaveResultsToplistCompletionLoop(WeaveResultsToplist *toplist)
Compute all missing 'extra' (non-toplist-ranking) statistics for all toplist entries.
Module which handles toplists of results.
REAL4(* WeaveResultsToplistItemGetRankStat)(const WeaveResultsToplistItem *item)
Function which returns the value of the statistic by which toplist items are ranked.
const REAL4 *(* WeaveResultsToplistRankingStats)(const WeaveSemiResults *semi_res)
Function which returns pointer to array of statistics by which toplist items are ranked.
void(* WeaveResultsToplistItemSetRankStat)(WeaveResultsToplistItem *item, const REAL4 value)
Function which sets the value of the statistic by which toplist items are ranked.
const char * name
Definition: SearchTiming.c:93
@ WEAVE_SIMULATE
Simulate search (implicitly with full memory allocation)
Definition: Weave.h:82
enum tagWeaveStatisticType WeaveStatisticType
Definition: Weave.h:61
int l
@ WEAVE_STATISTIC_MEAN2F_DET
Per detector average F-statistic.
@ WEAVE_STATISTIC_COH2F
Per segment multi-detector F-statistic.
@ WEAVE_STATISTIC_MAX2F
@ WEAVE_STATISTIC_NCOUNT_DET
Hough number count per detector.
@ WEAVE_STATISTIC_COH2F_DET
Per segment per-detector F-statistic.
@ WEAVE_STATISTIC_BtSGLtL
(transient-)line robust log10(B_tS/GLtL) statistic
@ WEAVE_STATISTIC_SUM2F_DET
Per detector sum F-statistic.
@ WEAVE_STATISTIC_MAX2F_DET
@ WEAVE_STATISTIC_MEAN2F
Multi-detector average (over segments) F-statistic.
@ WEAVE_STATISTIC_BSGLtL
(transient-)line robust log10(B_S/GLtL) statistic
@ WEAVE_STATISTIC_BSGL
Line-robust log10(B_S/GL) statistic.
@ WEAVE_STATISTIC_NCOUNT
Hough number count.
@ WEAVE_STATISTIC_SUM2F
Multi-detector sum (over segments) F-statistic.
#define WEAVE_STATISTIC_NAME(ws)
Names of all possible statistics.
const double u
const double r2
#define XLAL_FITS_TABLE_COLUMN_BEGIN(record_type)
Definition: FITSFileIO.h:239
#define XLAL_FITS_TABLE_COLUMN_ADD_NAMED(file, type, field, col_name)
Definition: FITSFileIO.h:246
#define XLAL_FITS_TABLE_COLUMN_ADD_PTR_ARRAY_NAMED(file, type, length, field, col_name)
Definition: FITSFileIO.h:264
struct tagFITSFile FITSFile
Representation of a FITS file.
Definition: FITSFileIO.h:54
unsigned char BOOLEAN
#define XLAL_NUM_ELEM(x)
uint64_t UINT8
double REAL8
#define XLAL_INIT_DECL(var,...)
uint32_t UINT4
float REAL4
int XLALHeapVisit(const LALHeap *h, LALHeapVisitFcn visit, void *visit_param)
int XLALHeapAdd(LALHeap *h, void **x)
int XLALHeapModify(LALHeap *h, LALHeapModifyFcn modify, void *modify_param)
int XLALHeapSize(const LALHeap *h)
const void * XLALHeapRoot(const LALHeap *h)
int XLALHeapIsFull(const LALHeap *h)
const void ** XLALHeapElements(const LALHeap *h)
void(* LALHeapDtorFcn)(void *x)
void XLALHeapDestroy(LALHeap *h)
LALHeap * XLALHeapCreate2(LALHeapDtorFcn dtor, int max_size, int min_or_max_heap, LALHeapCmpParamFcn cmp, void *cmp_param)
void * XLALCalloc(size_t m, size_t n)
void XLALFree(void *p)
#define LAL_UINT8_FORMAT
int XLALCompareREAL4Vectors(VectorComparison *result, const REAL4Vector *x, const REAL4Vector *y, const VectorComparison *tol)
Compare two REAL4 vectors using various different comparison metrics.
REAL4 XLALComputeBSGLtL(const REAL4 twoF, const REAL4 twoFX[PULSAR_MAX_DETECTORS], const REAL4 maxtwoFlX[PULSAR_MAX_DETECTORS], const BSGLSetup *setup)
Single-bin wrapper of XLALVectorComputeBSGLtL(), provided for backwards compatibility.
REAL4 XLALComputeBtSGLtL(const REAL4 maxtwoFl, const REAL4 twoFX[PULSAR_MAX_DETECTORS], const REAL4 maxtwoFXl[PULSAR_MAX_DETECTORS], const BSGLSetup *setup)
Single-bin wrapper of XLALVectorComputeBtSGLtL(), provided for backwards compatibility.
REAL4 XLALComputeBSGL(const REAL4 twoF, const REAL4 twoFX[PULSAR_MAX_DETECTORS], const BSGLSetup *setup)
Single-bin wrapper of XLALVectorComputeBSGL(), provided for backwards compatibility.
#define PULSAR_MAX_DETECTORS
maximal number of detectors we can handle (for static arrays of detector quantities)
#define PULSAR_MAX_SPINS
maximal number of spin-parameters (Freq + spindowns) we can handle
int XLALConvertPhysicalToSuperskyPoint(gsl_vector *out_rssky, const PulsarDopplerParams *in_phys, const SuperskyTransformData *rssky_transf)
Convert a point from physical to supersky coordinates.
REAL4Vector * XLALCreateREAL4Vector(UINT4 length)
void XLALDestroyUINT4Vector(UINT4Vector *vector)
UINT4Vector * XLALResizeUINT4Vector(UINT4Vector *vector, UINT4 length)
void XLALDestroyREAL4Vector(REAL4Vector *vector)
int XLALVectorFindScalarLessEqualREAL4(UINT4 *count, UINT4 *out, REAL4 scalar, const REAL4 *in, const UINT4 len)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_ERROR(...)
#define XLAL_CHECK(assertion,...)
#define XLAL_TRY(statement, errnum)
#define XLAL_CHECK_NULL(assertion,...)
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFUNC
XLAL_ETOL
XLAL_EERR
XLAL_EINVAL
XLAL_EFAILED
float data[BLOCKSIZE]
Definition: hwinject.c:360
list y
n
desc
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
PulsarSpins fkdot
Intrinsic spins: [Freq, f1dot, f2dot, ... ] where fkdot = d^kFreq/dt^k.
REAL8 Delta
Sky position: DEC (latitude) in equatorial coords and radians.
REAL8 Alpha
Sky position: RA (longitude) in equatorial coords and radians.
REAL4 * data
Struct holding the results of comparing two floating-point vectors (real-valued or complex),...
Definition: LFTandTSutils.h:64
Toplist of output results.
WeaveResultsToplistItemSetRankStat item_set_rank_stat_fcn
Function which sets the value of the statistic by which toplist items are ranked.
UINT4Vector * maybe_add_freq_idxs
Vector of indexes of toplist results which should be considered for addition.
WeaveResultsToplistItem * saved_item
Save a no-longer-used toplist item for re-use.
const char * stat_desc
Description of ranking statistic.
UINT8 serial
Serial number which is incremented for each item.
size_t nspins
Number of spindown parameters to output.
const char * stat_name
Name of ranking statistic.
WeaveResultsToplistItemGetRankStat item_get_rank_stat_fcn
Function which returns the value of the statistic by which toplist items are ranked.
WeaveStatisticsParams * statistics_params
Struct holding all parameters for which statistics to output and compute, when, and how NOTE: this is...
LALHeap * heap
Heap which ranks toplist items.
WeaveResultsToplistRankingStats rank_stats_fcn
Function which returns pointer to array of statistics by which toplist items are ranked.