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
splitSFTs.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2020, 2022 Karl Wette
3 * Copyright (C) 2008, 2010 Bernd Machenschalk, Bruce Allen
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with with program; see the file COPYING. If not, write to the
17 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
18 * MA 02110-1301 USA
19 */
20
21/**
22 * \file
23 * \ingroup lalpulsar_bin_SFTTools
24 * \author Bernd Machenschalk, Bruce Allen
25 *
26 * \brief This program reads in binary SFTs and writes out narrow-banded merged SFTs.
27 *
28 * Writen by Bernd Machenschalk for Einstein\@home 2008
29 *
30 * Revision splitSFTs.c,v 1.41 2008/10/29 16:54:13 was
31 * reviewed by the LSC CW Review Committee Thur Oct 30, 2008
32 *
33 * Suggested improvements:
34 * issue warning at ambiguous frequency values that are so close to the boundary of a bin
35 * that rounding might end up giving an unintended bin
36 * check for consistency of input SFTs (same timebase, ascending timestamps etc., see spec)
37 * and merged SFTs (check last header of a file we are appending to)
38 *
39 * Other possible improvements not suggested by the committee
40 * keep output files open (if there aren't too many)
41 * obscure a mystery factor in command-line record even if given with long option --factor
42 */
43
44#include "config.h"
45
46#include <math.h>
47#include <time.h>
48#ifdef HAVE_UNISTD_H
49#include <unistd.h>
50#endif
51#include <stdio.h>
52#include <stdlib.h>
53#include <string.h>
54#include <glob.h>
55#include <sys/stat.h>
56#include <lal/LALStdlib.h>
57#include <lal/LALHashTbl.h>
58#include <lal/Date.h>
59#include <lal/SFTfileIO.h>
60#include <lal/SFTReferenceLibrary.h>
61#include <lal/LALPulsarVCSInfo.h>
62
63#define FALSE 0
64#define TRUE (!FALSE)
65
66#define CMT_NONE 0
67#define CMT_OLD 1
68#define CMT_FULL 2
69
70// Compare two quantities, and return a sort order value if they are unequal
71#define COMPARE_BY( x, y ) do { if ( (x) < (y) ) return -1; if ( (x) > (y) ) return +1; } while(0)
72
73typedef struct {
74 int resource_units; /**< number of resource units available for consumption */
75 int resource_rate; /**< this many resource units become available in one second */
76 time_t last_checked; /**< time we last checked */
78
79/** throtteling settings */
84
85/** request a resurce. Function returns after waiting for throttle time */
86static void request_resource( UNIT_SOURCE *us, int units )
87{
88 time_t now;
89 int seconds;
90
91 /* negative rate indicates no throttling */
92 if ( us->resource_rate <= 0 ) {
93 return;
94 }
95
96 us->resource_units -= units;
97 while ( us->resource_units <= 0 ) {
98 time( &now );
99 seconds = now - us->last_checked;
100 /* guard against overflow and condor checkpointing */
101 if ( seconds < 0 ) {
102 seconds = 0;
103 }
104 if ( seconds > 10 ) {
105 seconds = 10;
106 }
107 if ( seconds == 0 ) {
108 sleep( 1 );
109 }
110 us->resource_units += seconds * us->resource_rate;
111 us->last_checked = now;
112 }
113}
114
115/* determine if the given filepath is an existing directory or not */
116static int is_directory( const char *fname )
117{
118 struct stat stat_buf;
119
120 if ( stat( fname, &stat_buf ) ) {
121 return 0;
122 }
123
124 if ( ! S_ISDIR( stat_buf.st_mode ) ) {
125 return 0;
126 } else {
127 return 1;
128 }
129
130}
131
132/* hash an SFTFilenameSpec */
133static UINT8 hash_SFTFilenameSpec( const void *x )
134{
135 const SFTFilenameSpec *rec = ( const SFTFilenameSpec * ) x;
136 UINT4 hval = 0;
137 XLALPearsonHash( &hval, sizeof( hval ), &rec->detector, sizeof( rec->detector ) );
138 XLALPearsonHash( &hval, sizeof( hval ), &rec->SFTtimebase, sizeof( rec->SFTtimebase ) );
139 XLALPearsonHash( &hval, sizeof( hval ), &rec->nbFirstBinFreq, sizeof( rec->nbFirstBinFreq ) );
140 XLALPearsonHash( &hval, sizeof( hval ), &rec->nbFirstBinRem, sizeof( rec->nbFirstBinRem ) );
141 XLALPearsonHash( &hval, sizeof( hval ), &rec->nbBinWidthFreq, sizeof( rec->nbBinWidthFreq ) );
142 XLALPearsonHash( &hval, sizeof( hval ), &rec->nbBinWidthRem, sizeof( rec->nbBinWidthRem ) );
143 return hval;
144}
145
146/* compare an SFTFilenameSpec */
147static int compare_SFTFilenameSpec( const void *x, const void *y )
148{
149 const SFTFilenameSpec *recx = ( const SFTFilenameSpec * ) x;
150 const SFTFilenameSpec *recy = ( const SFTFilenameSpec * ) y;
151 COMPARE_BY( recx->detector[0], recy->detector[0] );
152 COMPARE_BY( recx->detector[1], recy->detector[1] );
153 COMPARE_BY( recx->SFTtimebase, recy->SFTtimebase );
155 COMPARE_BY( recx->nbFirstBinRem, recy->nbFirstBinRem );
157 COMPARE_BY( recx->nbBinWidthRem, recy->nbBinWidthRem );
158 return 0;
159}
160
161static int move_to_next_SFT( FILE *fpin, int nsamples, int comment_length )
162{
163 int move;
164 move = sizeof( struct headertag2 ) + nsamples * 2 * sizeof( float ) + comment_length;
165 fseek( fpin, move, SEEK_CUR );
166 return XLAL_SUCCESS;
167}
168
169/** main program */
170int main( int argc, char **argv )
171{
172 int arg; /* current command-line argument */
173 unsigned int bin; /* current bin */
174 struct headertag2 hd, lasthd; /* header of input SFT */
175 FILE *fpin; /* currently open input filepointer */
176 FILE *fpout; /* currently open output filepointer */
177 char *oldcomment; /* comment of input SFT */
178 char *cmdline = NULL; /* records command-line to add it to comment */
179 char *comment = NULL; /* comment to be written into output SFT file */
180 int swap; /* do we need to swap bytes? */
181 float *data; /* SFT data */
182 char outdir0[] = ".";
183 char *outdir = ( char * )outdir0; /* output filename prefix */
184 double factor = 1.0; /* "mystery" factor */
185 int firstfile = TRUE; /* are we processing the first input SFT file? */
186 int allcomments = FALSE; /* write comment into _every_ SFT in the file */
187 int add_comment = CMT_FULL; /* add VCS ID and full command-line to every SFT file */
188 unsigned int startBin = 0, endBin = 0; /* start and end in bins */
189 unsigned int width = 0, overlap = 0; /* width and overlap in bins */
190 double fMin = -1.0, fMax = -1.0; /* start and end in Hz */
191 double fWidth = -1.0, fOverlap = -1; /* width and overlap in Hz */
192 unsigned int nactivesamples; /* number of bins to actually read in */
193 int validate = TRUE; /* validate the checksum of each input SFT before using it? */
194 LIGOTimeGPS *minStartTime = NULL; /* pointer for optional constraint */
195 LIGOTimeGPS *maxStartTime = NULL; /* pointer for optional constraint */
196 LIGOTimeGPS XLAL_INIT_DECL( userMinStartTime ); /* user-given constraint: earliest SFT timestamp to start using */
197 LIGOTimeGPS XLAL_INIT_DECL( userMaxStartTime ); /* user-given constraint: earliest SFT timestamp to no longer use */
198 int assumeSorted = 0; /* Are SFT input files chronologically sorted? */
199 int sfterrno = 0; /* SFT error number return from reference library */
200 LALHashTbl *nbsfts = NULL; /* hash table of existing narrow-band SFTs */
201
202 /* initialize throtteling */
207
208 /* help / usage message */
209 if ( ( argv[1] == NULL ) ||
210 ( strcmp( argv[1], "-h" ) == 0 ) ||
211 ( strcmp( argv[1], "--help" ) == 0 ) ) {
212 fprintf( stderr,
213 "%s -h\n"
214 "\n"
215 " Write this help message\n"
216 "\n"
217 "%s\n"
218 " [-a|--all-comments]\n"
219 " [-c|--add-comment 0|1|2]\n"
220 " [-v|--no-validation]\n"
221 " [-s|--start-bin <startbin>]\n"
222 " [-e|--end-bin <endbin (exclusively)>]\n"
223 " [-b|--width <sftbins>]\n"
224 " [-x|--overlap <overlap>]\n"
225 " [-fs|--start-frequency <startfrequency>]\n"
226 " [-fe|--end-frequency <endfrequency (exclusively)>]\n"
227 " [-fb|--frequency-bandwidth <frequencywidth>]\n"
228 " [-fx|--frequency-overlap <frequencyoverlap>]\n"
229 " [-ts|--minStartTime <minStartTime>]\n"
230 " [-te|--maxStartTime <maxStartTime (exclusively)>]\n"
231 " [-as|--assumeSorted 1|0]\n"
232 " [-m|--factor <factor>]\n"
233 " [-n|--output-directory <outputdirectory>]\n"
234 " [--] <inputfile> ...\n"
235 "\n"
236 " This program reads in binary SFTs and writes out narrow-banded\n"
237 " merged SFTs.\n"
238 "\n"
239 " The frequency bands of the ouput SFTs (first frequency bin of first output SFT,\n"
240 " last frequency bin of last output SFT, number of bins in each output SFT)\n"
241 " and a possible overlap of the output files can be specified\n"
242 " either in bins ('-s', '-e', '-b', -'x') or Hz ('-fs', '-fe', '-fb', '-fx')\n"
243 " (or mixed - if both are given, frequency values take precedence).\n"
244 "\n"
245 " A 'mystery factor' can be specified with '-m' option.\n"
246 "\n"
247 " '-v' skips validation of the CRC checksum of the input files\n"
248 "\n"
249 " The name of the output SFT follows the standard SFT filename convention, with\n"
250 " information specific to narrow-band SFTs included in the description:\n"
251 " <outdir>/<obs>-<nSFT>_<det>_<timebase>SFT_NBF<firstbin>W<binwidth>[_<misc>]-<start>-<span>.sft\n"
252 " where\n"
253 " <outdir> = output directory given by the '-n' option, default is '.'\n"
254 " <nSFT> = number of SFTs in the file\n"
255 " <obs> = 1-character observatory prefix, e.g. 'H', 'L', 'V'\n"
256 " <det> = 2-character detector prefix, e.g. 'H1', 'L1', 'V1'\n"
257 " <timebase> = time base of the SFTs, rounded to the nearest second\n"
258 " <firstbin> = first bin of the SFTs, expressed in the form\n"
259 " <firstbin> = <freq>Hz<remainder>, where\n"
260 " <freq> = <firstbin> / <timebase>, rounded down to nearest Hz\n"
261 " <remainder> = <firstbin> - <freq> * <timebase>, number of remaining bins\n"
262 " <binwidth> = bin width of the SFTs, expressed in the form\n"
263 " <binwidth> = <freq>Hz<remainder>, where\n"
264 " <freq> = <binwidth> / <timebase>, rounded down to nearest Hz\n"
265 " <remainder> = <binwidth> - <freq> * <timebase>, number of remaining bins\n"
266 " <misc> = misc description field, if any, from the input files\n"
267 " <start> = start time of the first SFT as a GPS timestamp\n"
268 " <span> = time spanned by the SFT in seconds\n"
269 " If an output file already exists, the program will append the new SFTs to them, making\n"
270 " it possible to construct the final narrow-band SFTs by running the program multiple\n"
271 " times with different input SFTs. The GPS timestamps of the input SFTs need to be in\n"
272 " ascending order to get valid merged SFT files.\n"
273 "\n"
274 " The '-c' options specifies how to deal with comments - 0 means no comment is written\n"
275 " at all, 1 means that the comment is taken unmodified from the input SFTs, 2 (default)\n"
276 " means that the program appends its RCS id and command-line to the comment.\n"
277 " By default a comment is written only to the first SFT of a merged SFT output 'block'\n"
278 " (i.e. call to this program). Adding the option '-a' to the command line specifies that\n"
279 " instead the comment is written into every SFT in the resulting file.\n"
280 "\n"
281 " The options 'minStartTime' and 'maxStartTime' constrain input SFTs to be read\n"
282 " according to their time stamps. Following the XLALCWGPSinRange() convention,\n"
283 " these are interpreted as half-open intervals:\n"
284 " \tminStartTime <= time stamp < maxStartTime.\n"
285 " After encountering a time stamp outside of the given range, the program jumps to\n"
286 " the next input file as specified in the input argument; i.e. time stamp constraints\n"
287 " are imposed file-wise. For a multi-SFT file, headers are checked to ensure proper\n"
288 " sorting. Mind that all timestamps refer to SFT **start** times.\n"
289 "\n"
290 " If 'assumeSorted' is set to 1, SFT input files will be assumed to be chronologically\n"
291 " sorted, which means the program will stop as soon as an SFT located after the\n"
292 " specified range is encountered.\n"
293 "\n"
294 " After all options (and an optional '--' separator), the input files are given, as many\n"
295 " as you wish (or the OS supports - using xargs should be simple with this command-line\n"
296 " syntax).\n"
297 "\n"
298 " The program adds its own VCS ID and command-line to the comment of the written SFTs,\n"
299 " a mystery factor should show up as 'xxx' there.\n",
300 argv[0], argv[0] );
301 exit( 0 );
302 }
303
304 /* record VCS ID and command-line for the comment */
305 XLAL_CHECK_MAIN( ( cmdline = ( char * )XLALMalloc( strlen( lalPulsarVCSIdentInfo.vcsId ) + strlen( lalPulsarVCSIdentInfo.vcsStatus ) + 2 ) ) != NULL, XLAL_ENOMEM, "out of memory allocating cmdline" );
306 strcpy( cmdline, lalPulsarVCSIdentInfo.vcsId );
307 strcat( cmdline, lalPulsarVCSIdentInfo.vcsStatus );
308 strcat( cmdline, "\n" );
309 for ( arg = 0; arg < argc; arg++ ) {
310 if ( strcmp( argv[arg], "-m" ) == 0 ) {
311 /* obscure the mystery factor */
312 XLAL_CHECK_MAIN( ( cmdline = ( char * )XLALRealloc( ( void * )cmdline, strlen( cmdline ) + 8 ) ) != NULL, XLAL_ENOMEM, "out of memory allocating cmdline" );
313 strcat( cmdline, "-m xxx " );
314 arg++;
315 } else {
316 XLAL_CHECK_MAIN( ( cmdline = ( char * )XLALRealloc( ( void * )cmdline, strlen( cmdline ) + strlen( argv[arg] ) + 2 ) ) != NULL, XLAL_ENOMEM, "out of memory allocating cmdline" );
317 strcat( cmdline, argv[arg] );
318 if ( arg == argc - 1 ) {
319 strcat( cmdline, "\n" );
320 } else {
321 strcat( cmdline, " " );
322 }
323 }
324 }
325
326 /* get parameters from command-line */
327 for ( arg = 1; arg < argc; arg++ ) {
328 if ( ( strcmp( argv[arg], "-c" ) == 0 ) ||
329 ( strcmp( argv[arg], "--add-comment" ) == 0 ) ) {
330 add_comment = atoi( argv[++arg] );
331 } else if ( ( strcmp( argv[arg], "-a" ) == 0 ) ||
332 ( strcmp( argv[arg], "--all-comments" ) == 0 ) ) {
333 allcomments = TRUE;
334 } else if ( ( strcmp( argv[arg], "-s" ) == 0 ) ||
335 ( strcmp( argv[arg], "--start-bin" ) == 0 ) ) {
336 startBin = atoi( argv[++arg] );
337 } else if ( ( strcmp( argv[arg], "-e" ) == 0 ) ||
338 ( strcmp( argv[arg], "--end-bin" ) == 0 ) ) {
339 endBin = atoi( argv[++arg] );
340 } else if ( ( strcmp( argv[arg], "-b" ) == 0 ) ||
341 ( strcmp( argv[arg], "--width" ) == 0 ) ) {
342 width = atoi( argv[++arg] );
343 } else if ( ( strcmp( argv[arg], "-x" ) == 0 ) ||
344 ( strcmp( argv[arg], "--overlap" ) == 0 ) ) {
345 overlap = atoi( argv[++arg] );
346 } else if ( ( strcmp( argv[arg], "-fs" ) == 0 ) ||
347 ( strcmp( argv[arg], "--start-frequency" ) == 0 ) ) {
348 fMin = atof( argv[++arg] );
349 } else if ( ( strcmp( argv[arg], "-fe" ) == 0 ) ||
350 ( strcmp( argv[arg], "--end-frequency" ) == 0 ) ) {
351 fMax = atof( argv[++arg] );
352 } else if ( ( strcmp( argv[arg], "-fb" ) == 0 ) ||
353 ( strcmp( argv[arg], "--frequency-bandwidth" ) == 0 ) ) {
354 fWidth = atof( argv[++arg] );
355 } else if ( ( strcmp( argv[arg], "-fx" ) == 0 ) ||
356 ( strcmp( argv[arg], "--frequency-overlap" ) == 0 ) ) {
357 fOverlap = atof( argv[++arg] );
358 } else if ( ( strcmp( argv[arg], "-ts" ) == 0 ) ||
359 ( strcmp( argv[arg], "--minStartTime" ) == 0 ) ) {
360 XLAL_CHECK_MAIN( XLALGPSSetREAL8( &userMinStartTime, ( REAL8 )atof( argv[++arg] ) ) != NULL, XLAL_EDOM );
361 minStartTime = &userMinStartTime;
362 } else if ( ( strcmp( argv[arg], "-te" ) == 0 ) ||
363 ( strcmp( argv[arg], "--maxStartTime" ) == 0 ) ) {
364 XLAL_CHECK_MAIN( XLALGPSSetREAL8( &userMaxStartTime, ( REAL8 )atof( argv[++arg] ) ) != NULL, XLAL_EDOM );
365 maxStartTime = &userMaxStartTime;
366 } else if ( ( strcmp( argv[arg], "-as" ) == 0 ) ||
367 ( strcmp( argv[arg], "--assumeSorted" ) == 0 ) ) {
368 assumeSorted = atoi( argv[++arg] );
369 } else if ( ( strcmp( argv[arg], "-m" ) == 0 ) ||
370 ( strcmp( argv[arg], "--factor" ) == 0 ) ) {
371 factor = atof( argv[++arg] );
372 } else if ( ( strcmp( argv[arg], "-v" ) == 0 ) ||
373 ( strcmp( argv[arg], "--no-validation" ) == 0 ) ) {
374 validate = FALSE;
375 } else if ( ( strcmp( argv[arg], "-n" ) == 0 ) ||
376 ( strcmp( argv[arg], "--output-directory" ) == 0 ) ) {
377 outdir = argv[++arg];
378 } else if ( ( strcmp( argv[arg], "-rb" ) == 0 ) ||
379 ( strcmp( argv[arg], "--read-bandwidth" ) == 0 ) ) {
380 read_bandwidth.resource_rate = atoi( argv[++arg] );
381 } else if ( ( strcmp( argv[arg], "-ror" ) == 0 ) ||
382 ( strcmp( argv[arg], "--read-open-rate" ) == 0 ) ) {
383 read_open_rate.resource_rate = atoi( argv[++arg] );
384 } else if ( ( strcmp( argv[arg], "-wb" ) == 0 ) ||
385 ( strcmp( argv[arg], "--write-bandwidth" ) == 0 ) ) {
386 write_bandwidth.resource_rate = atoi( argv[++arg] );
387 } else if ( ( strcmp( argv[arg], "-wor" ) == 0 ) ||
388 ( strcmp( argv[arg], "--write-open-rate" ) == 0 ) ) {
389 write_open_rate.resource_rate = atoi( argv[++arg] );
390 } else if ( strcmp( argv[arg], "--" ) == 0 ) {
391 ++arg;
392 break;
393 } else if ( strncmp( argv[arg], "-", 1 ) == 0 ) {
394 fprintf( stderr, "unknown option '%s', try '-h' for help\n", argv[arg] );
395 exit( -1 );
396 } else {
397 break;
398 }
399 }
400
401 /* build an output string for the user timestamps constraint */
402 char *constraint_str;
403 const size_t constraint_str_len = 32; /* just needs to be big enough for two INT4 and a comma */
404 XLAL_CHECK_MAIN( ( constraint_str = ( char * )XLALMalloc( constraint_str_len ) ) != NULL,
405 XLAL_ENOMEM, "out of memory allocating constraint_str" );
406 if ( minStartTime ) {
407 XLAL_CHECK_MAIN( snprintf( constraint_str, constraint_str_len, "%d,", minStartTime->gpsSeconds ) < ( int )constraint_str_len,
408 XLAL_ESYS, "output of snprintf() was truncated" );
409 } else {
410 XLAL_CHECK_MAIN( snprintf( constraint_str, constraint_str_len, "," ) < ( int )constraint_str_len,
411 XLAL_ESYS, "output of snprintf() was truncated" );
412 }
413 if ( maxStartTime ) {
414 XLAL_CHECK_MAIN( snprintf( constraint_str + strlen( constraint_str ), constraint_str_len - strlen( constraint_str ), "%d", maxStartTime->gpsSeconds ) < ( int )( constraint_str_len - strlen( constraint_str ) ),
415 XLAL_ESYS, "output of snprintf() was truncated" );
416 }
417
418 /* check if there was an input-file option given at all */
419 XLAL_CHECK_MAIN( argv[arg] != NULL, XLAL_EINVAL, "no input files specified" );
420
421 /* check output directory exists */
422 XLAL_CHECK_MAIN( is_directory( outdir ), XLAL_ESYS, "output directory does not exist" );
423
424 /* find existing narrow-band SFTs */
426 XLAL_CHECK_MAIN( nbsfts != NULL, XLAL_EFUNC, "XLALHashTblCreate() failed" );
427 {
428 char *pattern = XLALStringAppendFmt( NULL, "%s/*.sft", outdir );
429 XLAL_CHECK( pattern != NULL, XLAL_ENOMEM );
430 glob_t globbuf;
431 glob( pattern, GLOB_MARK | GLOB_NOSORT, NULL, &globbuf );
432 for ( size_t i = 0; i < globbuf.gl_pathc; ++i ) {
433
434 /* parse SFT filename */
435 SFTFilenameSpec nbSFTspec;
436 XLAL_CHECK_MAIN( XLALParseSFTFilenameIntoSpec( &nbSFTspec, globbuf.gl_pathv[i] ) == XLAL_SUCCESS, XLAL_EFUNC );
437
438 /* filter the narrow-band SFTs */
439 if ( nbSFTspec.nbFirstBinFreq == 0 ) {
440 continue;
441 }
442
443 /* create new SFT record */
444 SFTFilenameSpec *rec = NULL;
445 XLAL_CHECK_MAIN( ( rec = XLALMalloc( sizeof( *rec ) ) ) != NULL, XLAL_ENOMEM, "out of memory allocating SFT record" );
446 *rec = nbSFTspec;
447
448 /* add SFT record */
449 XLALPrintInfo( "Existing SFT: det=%c%c timebase=%u firstbin=(%u,%u) binwidth=(%u,%u) filename=%s\n",
450 rec->detector[0], rec->detector[1], rec->SFTtimebase,
452 globbuf.gl_pathv[i] );
453 XLAL_CHECK_MAIN( XLALHashTblAdd( nbsfts, rec ) == XLAL_SUCCESS, XLAL_EFUNC, "XLALHashTblAdd() failed" );
454
455 }
456 globfree( &globbuf );
457 XLALFree( pattern );
458 }
459
460 int stopInputCauseSorted = 0;
461 /* loop over all input SFT files */
462 for ( ; ( arg < argc ) && !stopInputCauseSorted; arg++ ) {
463
464 /* parse SFT filename */
465 SFTFilenameSpec inputSFTspec;
466 XLAL_CHECK_MAIN( XLALParseSFTFilenameIntoSpec( &inputSFTspec, argv[arg] ) == XLAL_SUCCESS, XLAL_EFUNC );
467
468 /* open input SFT file */
470 XLAL_CHECK_MAIN( ( fpin = fopen( argv[arg], "r" ) ) != NULL, XLAL_EIO, "could not open SFT file for reading" );
471
472 /* loop until end of input SFT file, when ReadSFTHeader() will return SFTENONE */
473 int firstread = TRUE;
474 int nSFT_this_file = 0;
475 while ( 1 ) {
476
477 /* read header, break if ReadSFTHeader() returns SFTENONE
478 (unless this is the first read, i.e. fail is SFT file is completely empty) */
480 sfterrno = ReadSFTHeader( fpin, &hd, &oldcomment, &swap, validate );
481 if ( !firstread ) {
482 if ( sfterrno == SFTENONE ) {
483 break;
484 }
485 }
486 XLAL_CHECK_MAIN( sfterrno == 0, XLAL_EIO, "could not read SFT header: %s", SFTErrorMessage( sfterrno ) );
487
488 /* Check that various bits of header information are consistent.
489 * This includes a check for monotonic increasing timestamps.
490 */
491 XLAL_CHECK_MAIN( firstread || !( sfterrno = CheckSFTHeaderConsistency( &lasthd, &hd ) ), XLAL_EIO, "Inconsistent SFT headers: %s", SFTErrorMessage( sfterrno ) );
492 lasthd = hd; /* keep copy of header for comparison the next time */
493
494 /* only use SFTs within specified ranges */
495 LIGOTimeGPS startTimeGPS = {hd.gps_sec, 0};
496 int inRange = XLALCWGPSinRange( startTimeGPS, minStartTime, maxStartTime );
497 XLAL_CHECK_MAIN( !( firstread && ( inRange == 1 ) ), XLAL_EIO, "First timestamp %d in file was after user constraint [%s)!", hd.gps_sec, constraint_str );
498 firstread = FALSE;
499 if ( inRange == -1 ) { /* input timestamp too early, skip */
501 continue;
502 }
503 if ( inRange == 1 ) { /* input timestamp too late, stop with this file */
504 stopInputCauseSorted = assumeSorted;
505 break;
506 }
507 nSFT_this_file += 1;
508
509 /* calculate bins from frequency parameters if they were given */
510 /* deltaF = 1.0 / tbase; bins = freq / deltaF => bins = freq * tbase */
511 {
512 const double deltaF = 1.0 / hd.tbase;
513 if ( fMin >= 0.0 ) {
515 }
516 if ( fMax >= 0.0 ) {
517 endBin = XLALRoundFrequencyUpToSFTBin( fMax, deltaF ) - 1;
518 }
519 if ( fWidth >= 0.0 ) {
520 width = XLALRoundFrequencyUpToSFTBin( fWidth, deltaF );
521 }
522 if ( fOverlap >= 0.0 ) {
523 overlap = XLALRoundFrequencyUpToSFTBin( fOverlap, deltaF );
524 }
525 }
526
527 /* allocate space for SFT data */
528 XLAL_CHECK_MAIN( ( data = ( float * )XLALCalloc( hd.nsamples, 2 * sizeof( float ) ) ) != NULL, XLAL_ENOMEM, "out of memory allocating data" );
529
530 /* error if desired start bin < hd.firstfreqindex */
531 if ( ( int )startBin < hd.firstfreqindex ) {
532 fprintf( stderr,
533 "ERROR: start bin (%d) is smaller than first bin in input SFT (%d)\n",
534 startBin, hd.firstfreqindex );
535 exit( 9 );
536 }
537
538 /* error if desired end bin > hd.firstfreqindex + hd.nsamples - 1 */
539 if ( startBin + width > endBin + 1 ) {
540 fprintf( stderr,
541 "ERROR: end bin (%d) is larger than last bin in input SFT (%d)\n",
542 endBin, hd.firstfreqindex + hd.nsamples - 1 );
543 exit( 10 );
544 }
545
546 /* error if overlap is larger than the width */
547 if ( overlap >= width ) {
548 fprintf( stderr,
549 "ERROR: overlap (%d) is not smaller than the width (%d)\n",
550 overlap, width );
551 exit( 11 );
552 }
553
554 /* construct comment for output SFTs */
555 if ( add_comment > CMT_OLD ) {
556
557 /* allocate space for new comment */
558 XLAL_CHECK_MAIN( ( comment = ( char * )XLALMalloc( hd.comment_length + strlen( cmdline ) + 1 ) ) != NULL, XLAL_ENOMEM, "out of memory allocating comment" );
559
560 /* append the commandline of this program to the old comment */
561 if ( oldcomment ) {
562 strcpy( comment, oldcomment );
563 } else {
564 *comment = '\0';
565 }
566 strcat( comment, cmdline );
567
568 } else if ( add_comment == CMT_OLD ) {
569
570 /* only copied existing comment, no additional space needed */
571 comment = oldcomment;
572
573 } /* else (add_comment == CMT_NONE) and (comment == NULL) i.e. no comment at all */
574
575 /* get the detector name from SFT header */
576 const char detector[] = { hd.detector[0], hd.detector[1], 0 };
577
578 /* calculate number of bins to actually read (from width + overlap) */
579 /* add width-overlap samples as lon as they are < the total number og bins to write */
580 for ( nactivesamples = 0; nactivesamples < endBin - startBin; nactivesamples += width - overlap );
581 /* add the last overlap */
582 nactivesamples += overlap + 1;
583 /* if this number is larger than the bins in the input sft, just use the latter */
584 if ( nactivesamples > hd.nsamples + hd.firstfreqindex - startBin ) {
585 nactivesamples = hd.nsamples + hd.firstfreqindex - startBin;
586 }
587
588 /* read in SFT bins */
589 request_resource( &read_bandwidth, nactivesamples * 8 );
590 sfterrno = ReadSFTData( fpin, data, startBin, nactivesamples, NULL, NULL );
591 XLAL_CHECK_MAIN( sfterrno == 0, XLAL_EIO, "could not read SFT data: %s", SFTErrorMessage( sfterrno ) );
592
593 /* apply mystery factor and possibly normalization factor */
594 for ( bin = 0; bin < 2 * nactivesamples; bin++ ) {
595 data[bin] *= factor;
596 }
597
598 /* loop over start bins for output SFTs */
599 for ( bin = startBin; bin < endBin; bin += width - overlap ) {
600 /* determine the number of bins actually to write from the desired 'width',
601 given that the remaining number of bin may be odd (especially from overlapping)
602 and the bins to write need to be present in the input sft
603 */
604 int last_input_bin = hd.firstfreqindex + hd.nsamples - 1;
605 int last_output_bin = bin + width - 1;
606 int max_input_width = last_output_bin <= last_input_bin ? width : width - ( last_output_bin - last_input_bin );
607 int max_output_width = endBin - bin + 1;
608 int this_width = max_input_width < max_output_width ? max_input_width : max_output_width;
609
610 int timebase = ( int ) round( hd.tbase );
611 int outfreq = ( int )floor( bin / timebase );
612 int outfreqbin = bin - outfreq * timebase;
613 int outwidth = ( int )floor( this_width / timebase );
614 int outwidthbin = this_width - outwidth * timebase;
615
616 /* check if this narrow-band SFT exists */
618 .detector = { detector[0], detector[1], 0 },
619 .SFTtimebase = timebase,
620 .nbFirstBinFreq = outfreq,
621 .nbFirstBinRem = outfreqbin,
622 .nbBinWidthFreq = outwidth,
623 .nbBinWidthRem = outwidthbin,
624 };
625 XLALPrintInfo( "Looking for SFT:det=%c%c timebase=%u firstbin=(%u,%u) binwidth=(%u,%u)\n",
626 key.detector[0], key.detector[1], key.SFTtimebase,
627 key.nbFirstBinFreq, key.nbFirstBinRem, key.nbBinWidthFreq, key.nbBinWidthRem );
628 SFTFilenameSpec *rec = NULL;
629 XLAL_CHECK_MAIN( XLALHashTblExtract( nbsfts, &key, ( void ** ) &rec ) == XLAL_SUCCESS, XLAL_EFUNC, "XLALHashTblExtract() failed" );
630 char *existingSFTfilename = NULL;
631 if ( rec == NULL ) {
632
633 /* create new SFT record, starting from input SFT filename */
634 XLAL_CHECK_MAIN( ( rec = XLALCalloc( 1, sizeof( *rec ) ) ) != NULL, XLAL_ENOMEM, "out of memory allocating SFT record" );
635 *rec = inputSFTspec;
636
637 /* change to output directory */
638 XLAL_CHECK_MAIN( XLALFillSFTFilenameSpecStrings( rec, outdir, NULL, NULL, NULL, NULL, NULL, NULL ) == XLAL_SUCCESS, XLAL_EFUNC );
639
640 /* zero number of SFTs & SFT span (to be updated) */
641 rec->numSFTs = 0;
642 rec->SFTspan = 0;
643
644 /* set starting GPS time from header (respecting -ts, -te) */
645 rec->gpsStart = hd.gps_sec;
646
647 /* add narrow-band fields */
648 rec->nbFirstBinFreq = outfreq;
649 rec->nbFirstBinRem = outfreqbin;
650 rec->nbBinWidthFreq = outwidth;
651 rec->nbBinWidthRem = outwidthbin;
652
653 } else {
654
655 /* reconstruct filename of existing SFT */
656 existingSFTfilename = XLALBuildSFTFilenameFromSpec( rec );
657 XLAL_CHECK( existingSFTfilename != NULL, XLAL_EFUNC );
658
659 /* check if added record would break increasing timestamps */
660 int rec_last_gpsStart = rec->gpsStart + rec->SFTspan - rec->SFTtimebase;
661 XLAL_CHECK_MAIN( hd.gps_sec >= rec_last_gpsStart, XLAL_EDOM,
662 "New SFT timestamp %d is before startTime+span-timebase=%u+%u-%u=%u from existing file %s. Appending would yield an invalid merged SFT file.",
663 hd.gps_sec, rec->gpsStart, rec->SFTspan, rec->SFTtimebase, rec_last_gpsStart, existingSFTfilename );
664
665 }
666
667 /* update number of SFTs and SFT timespan */
668 rec->numSFTs += 1;
669 rec->SFTspan = hd.gps_sec - rec->gpsStart + rec->SFTtimebase;
670 if ( hd.gps_nsec > 0 ) {
671 rec->SFTspan += 1;
672 }
673
674 /* construct filename for this output SFT */
675 char *newSFTfilename = XLALBuildSFTFilenameFromSpec( rec );
676 XLAL_CHECK( newSFTfilename != NULL, XLAL_EFUNC );
677
678 /* update SFT filename */
679 if ( existingSFTfilename != NULL && strcmp( existingSFTfilename, newSFTfilename ) != 0 ) {
680 XLALPrintInfo( "Renaming SFT '%s' to '%s'\n", existingSFTfilename, newSFTfilename );
681 rename( existingSFTfilename, newSFTfilename );
682 }
683
684 /* store new SFT record */
685 XLALPrintInfo( "New/updated SFT: det=%c%c timebase=%u firstbin=(%u,%u) binwidth=(%u,%u) filename=%s\n",
686 rec->detector[0], rec->detector[1], rec->SFTtimebase,
688 newSFTfilename );
689 XLAL_CHECK_MAIN( XLALHashTblAdd( nbsfts, rec ) == XLAL_SUCCESS, XLAL_EFUNC, "XLALHashTblAdd() failed" );
690
691 /* append this SFT */
693 XLAL_CHECK_MAIN( ( fpout = fopen( newSFTfilename, "a" ) ) != NULL, XLAL_EIO, "could not open SFT for writing" );
694
695 /* write the data */
696 /* write the comment only to the first SFT of a "block", i.e. of a call of this program */
697 request_resource( &write_bandwidth, 40 + this_width * 8 );
698 sfterrno = WriteSFT( fpout, hd.gps_sec, hd.gps_nsec, hd.tbase, bin, this_width, detector, hd.windowspec, ( firstfile || allcomments ) ? comment : NULL, data + 2 * ( bin - startBin ) );
699 XLAL_CHECK_MAIN( sfterrno == 0, XLAL_EIO, "could not write SFT data: %s", SFTErrorMessage( sfterrno ) );
700
701 /* close output SFT file */
702 fclose( fpout );
703
704 XLALFree( existingSFTfilename );
705 XLALFree( newSFTfilename );
706
707 } /* loop over output SFTs */
708
709 /* cleanup */
710 if ( add_comment > CMT_OLD ) {
711 XLALFree( comment );
712 }
713 XLALFree( data );
714
715 /* next file is not the first file anymore */
716 firstfile = FALSE;
717
718 /* Move forward to next SFT in merged file */
720
721 } /* end loop over SFTs in this file */
722
723 XLAL_CHECK_MAIN( nSFT_this_file > 0, XLAL_EIO, "Found no matching SFTs in file." );
724
725 /* close the input SFT file */
726 fclose( fpin );
727
728 } /* loop over input SFT files */
729
730 /* cleanup */
731 XLALFree( constraint_str );
732 if ( add_comment > CMT_OLD ) {
733 XLALFree( cmdline );
734 }
735 XLALHashTblDestroy( nbsfts );
737
738 return ( 0 );
739}
void LALCheckMemoryLeaks(void)
const LALVCSInfo lalPulsarVCSIdentInfo
Identable VCS and build information for LALPulsar.
const char * comment
Definition: SearchTiming.c:94
#define fprintf
uint64_t UINT8
double REAL8
#define XLAL_INIT_DECL(var,...)
uint32_t UINT4
int XLALPearsonHash(void *hval, const size_t hval_len, const void *data, const size_t data_len)
int XLALHashTblAdd(LALHashTbl *ht, void *x)
int XLALHashTblExtract(LALHashTbl *ht, const void *x, void **y)
LALHashTbl * XLALHashTblCreate(LALHashTblDtorFcn dtor, LALHashTblHashFcn hash, LALHashTblCmpFcn cmp)
void XLALHashTblDestroy(LALHashTbl *ht)
void * XLALMalloc(size_t n)
void * XLALCalloc(size_t m, size_t n)
void * XLALRealloc(void *p, size_t n)
void XLALFree(void *p)
char * XLALStringAppendFmt(char *s, const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(2
const char * SFTErrorMessage(int errorcode)
#define SFTENONE
int CheckSFTHeaderConsistency(struct headertag2 *headerone, struct headertag2 *headertwo)
int ReadSFTHeader(FILE *fp, struct headertag2 *info, char **comment, int *swapendian, int validate)
int WriteSFT(FILE *fp, int gps_sec, int gps_nsec, double tbase, int firstfreqindex, int nsamples, const char *detector, unsigned short windowspec, const char *comment, float *data)
int ReadSFTData(FILE *fp, float *data, int firstbin, int nsamples, char **comment, struct headertag2 *info)
char * XLALBuildSFTFilenameFromSpec(const SFTFilenameSpec *spec)
Build an SFT file name from the given specification.
Definition: SFTnaming.c:302
UINT4 XLALRoundFrequencyUpToSFTBin(const REAL8 freq, const REAL8 df)
Round a REAL8 frequency up to the nearest integer SFT bin number.
Definition: SFTtypes.c:82
UINT4 XLALRoundFrequencyDownToSFTBin(const REAL8 freq, const REAL8 df)
Round a REAL8 frequency down to the nearest integer SFT bin number.
Definition: SFTtypes.c:70
int XLALCWGPSinRange(const LIGOTimeGPS gps, const LIGOTimeGPS *minGPS, const LIGOTimeGPS *maxGPS)
Defines the official CW convention for whether a GPS time is 'within' a given range,...
Definition: SFTtypes.c:52
int XLALParseSFTFilenameIntoSpec(SFTFilenameSpec *spec, const char *SFTpath)
Parse a SFT file path and return its specification.
Definition: SFTnaming.c:463
int XLALFillSFTFilenameSpecStrings(SFTFilenameSpec *spec, const CHAR *path, const CHAR *extn, const CHAR *detector, const CHAR *window_type, const CHAR *privMisc, const CHAR *pubObsKind, const CHAR *pubChannel)
Convenience function for filling out the string fields in a SFTFilenameSpec.
Definition: SFTnaming.c:257
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_MAIN(assertion,...)
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EDOM
XLAL_EIO
XLAL_ESYS
XLAL_EINVAL
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
float data[BLOCKSIZE]
Definition: hwinject.c:360
list y
key
def now
int deltaF
static UINT8 hash_SFTFilenameSpec(const void *x)
Definition: splitSFTs.c:133
int main(int argc, char **argv)
main program
Definition: splitSFTs.c:170
static int is_directory(const char *fname)
Definition: splitSFTs.c:116
#define COMPARE_BY(x, y)
Definition: splitSFTs.c:71
UNIT_SOURCE read_bandwidth
throtteling settings
Definition: splitSFTs.c:80
static void request_resource(UNIT_SOURCE *us, int units)
request a resurce.
Definition: splitSFTs.c:86
UNIT_SOURCE write_open_rate
Definition: splitSFTs.c:83
#define TRUE
Definition: splitSFTs.c:64
#define FALSE
Definition: splitSFTs.c:63
static int compare_SFTFilenameSpec(const void *x, const void *y)
Definition: splitSFTs.c:147
#define CMT_OLD
Definition: splitSFTs.c:67
#define CMT_FULL
Definition: splitSFTs.c:68
UNIT_SOURCE write_bandwidth
Definition: splitSFTs.c:82
static int move_to_next_SFT(FILE *fpin, int nsamples, int comment_length)
Definition: splitSFTs.c:161
UNIT_SOURCE read_open_rate
Definition: splitSFTs.c:81
const char *const vcsStatus
const char *const vcsId
Structure specifying an SFT file name, following the convention in .
Definition: SFTfileIO.h:266
UINT4 SFTtimebase
Timebase in seconds of the SFT; set by XLALWriteSFT[Vector]2StandardFile()
Definition: SFTfileIO.h:271
UINT4 SFTspan
Total time interval in seconds covered by SFT file; set by XLALWriteSFT[Vector]2StandardFile()
Definition: SFTfileIO.h:275
UINT4 nbBinWidthRem
For narrow-band SFTs: remainder of division of SFT bandwidth by SFT time base.
Definition: SFTfileIO.h:284
UINT4 gpsStart
GPS time in seconds at the beginning of the first SFT in the file; set by XLALWriteSFT[Vector]2Standa...
Definition: SFTfileIO.h:274
UINT4 nbBinWidthFreq
For narrow-band SFTs: SFT bandwidth divided by SFT time base, rounded down.
Definition: SFTfileIO.h:283
CHAR detector[3]
2-character detector prefix (e.g.
Definition: SFTfileIO.h:270
UINT4 nbFirstBinRem
For narrow-band SFTs: remainder of division of SFT first bin frequency by SFT time base.
Definition: SFTfileIO.h:282
UINT4 numSFTs
Number of SFTs in the file; set by XLALWriteSFT[Vector]2StandardFile()
Definition: SFTfileIO.h:269
UINT4 nbFirstBinFreq
For narrow-band SFTs: SFT first bin frequency divided by SFT time base, rounded down.
Definition: SFTfileIO.h:281
int resource_rate
this many resource units become available in one second
Definition: splitSFTs.c:75
int resource_units
number of resource units available for consumption
Definition: splitSFTs.c:74
time_t last_checked
time we last checked
Definition: splitSFTs.c:76
unsigned short windowspec