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
compareSFTs.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2007 Reinhard Prix
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_SFTTools
23 * \author R. Prix
24 */
25
26/* ---------- includes ---------- */
27#include "config.h"
28
29#include <lal/Date.h>
30#include <lal/UserInput.h>
31#include <lal/SFTfileIO.h>
32#include <lal/LALPulsarVCSInfo.h>
33
34/* User variables */
35typedef struct {
36 CHAR *sftBname1;
37 CHAR *sftBname2;
38 INT4 debug;
40 BOOLEAN quiet;
41 REAL8 relErrorMax;
43
44/* local prototypes */
45int initUserVars( UserInput_t *uvar );
46REAL4 getMaxErrSFT( const SFTtype *sft1, const SFTtype *sft2 );
47REAL4 getMaxErrSFTVector( const SFTVector *sftvect1, const SFTVector *sftvect2 );
48REAL4 scalarProductSFT( const SFTtype *sft1, const SFTtype *sft2 );
49REAL4 scalarProductSFTVector( const SFTVector *sftvect1, const SFTVector *sftvect2 );
50SFTVector *subtractSFTVectors( const SFTVector *sftvect1, const SFTVector *sftvect2 );
51
52/*----------------------------------------------------------------------
53 * main function
54 *----------------------------------------------------------------------*/
55int
56main( int argc, char *argv[] )
57{
58 SFTVector *diffs = NULL;
59 REAL8 maxd = 0;
60
61 /* register all user-variables */
64
65 /* read cmdline & cfgfile */
66 BOOLEAN should_exit = 0;
68 if ( should_exit ) {
69 exit( 1 );
70 }
71
72 SFTCatalog *catalog;
73 XLAL_CHECK_MAIN( ( catalog = XLALSFTdataFind( uvar.sftBname1, NULL ) ) != NULL, XLAL_EFUNC );
74 SFTVector *SFTs1;
75 XLAL_CHECK_MAIN( ( SFTs1 = XLALLoadSFTs( catalog, -1, -1 ) ) != NULL, XLAL_EFUNC );
76 XLALDestroySFTCatalog( catalog );
77
78 XLAL_CHECK_MAIN( ( catalog = XLALSFTdataFind( uvar.sftBname2, NULL ) ) != NULL, XLAL_EFUNC );
79 SFTVector *SFTs2;
80 XLAL_CHECK_MAIN( ( SFTs2 = XLALLoadSFTs( catalog, -1, -1 ) ) != NULL, XLAL_EFUNC );
81 XLALDestroySFTCatalog( catalog );
82
83 /* ---------- do some sanity checks of consistency of SFTs ----------*/
84 XLAL_CHECK_MAIN( SFTs1->length == SFTs2->length, XLAL_EINVAL, "Number of SFTs differ for SFTbname1 and SFTbname2!\n" );
85
86 for ( UINT4 i = 0; i < SFTs1->length; i++ ) {
87 SFTtype *sft1 = &( SFTs1->data[i] );
88 SFTtype *sft2 = &( SFTs2->data[i] );
89
90 XLAL_CHECK_MAIN( strcmp( sft1->name, sft2->name ) == 0, XLAL_EINVAL, "\nERROR SFT %d: detector-prefix differ! '%s' != '%s'\n", i, sft1->name, sft2->name );
91
92 XLAL_CHECK_MAIN( sft1->data->length == sft2->data->length, XLAL_EINVAL, "\nERROR SFT %d: lengths differ! %d != %d\n", i, sft1->data->length, sft2->data->length );
93
94 REAL8 Tdiff = XLALGPSDiff( &( sft1->epoch ), &( sft2->epoch ) );
95 CHAR buf1[32], buf2[32];;
96 XLAL_CHECK_MAIN( Tdiff == 0.0, XLAL_EINVAL, "SFT %d: epochs differ: %s vs %s\n", i, XLALGPSToStr( buf1, &sft1->epoch ), XLALGPSToStr( buf2, &sft2->epoch ) );
97 XLAL_CHECK_MAIN( sft1->f0 == sft2->f0, XLAL_EINVAL, "ERROR SFT %d: fmin differ: %fHz vs %fHz\n", i, sft1->f0, sft2->f0 );
98 XLAL_CHECK_MAIN( sft1->deltaF == sft2->deltaF, XLAL_EINVAL, "ERROR SFT %d: deltaF differs: %fHz vs %fHz\n", i, sft1->deltaF, sft2->deltaF );
99 } /* for i < numSFTs */
100
101 /*---------- now do some actual comparisons ----------*/
102 XLAL_CHECK_MAIN( ( diffs = subtractSFTVectors( SFTs1, SFTs2 ) ) != NULL, XLAL_EFUNC );
103
104 if ( uvar.verbose ) {
105 for ( UINT4 i = 0; i < SFTs1->length; i++ ) {
106 SFTtype *sft1 = &( SFTs1->data[i] );
107 SFTtype *sft2 = &( SFTs2->data[i] );
108
109 XLALPrintInfo( "i=%02d: ", i );
110 REAL4 norm1 = scalarProductSFT( sft1, sft1 );
111 norm1 = sqrt( norm1 );
112 REAL4 norm2 = scalarProductSFT( sft2, sft2 );
113 norm2 = sqrt( norm2 );
114 REAL4 scalar = scalarProductSFT( sft1, sft2 );
115
116 REAL4 normdiff = scalarProductSFT( &( diffs->data[i] ), &( diffs->data[i] ) );
117
118 REAL4 d1 = ( norm1 - norm2 ) / norm1;
119 REAL4 d2 = 1.0 - scalar / ( norm1 * norm2 );
120 REAL4 d3 = normdiff / ( norm1 * norm1 + norm2 * norm2 );
121 REAL4 d4 = getMaxErrSFT( sft1, sft2 );
122 XLALPrintInfo( "SFT #%d: (|x|-|y|)/|x|=%10.3e, 1-x.y/(|x||y|)=%10.3e, |x-y|^2/(|x|^2+|y|^2))=%10.3e, maxErrSFT=%10.3e\n", i, d1, d2, d3, d4 );
123 } /* for i < SFTs->length */
124 } /* if verbose */
125
126 /* ---------- COMBINED measures ---------- */
127 {
128 REAL4 ret;
129 ret = scalarProductSFTVector( SFTs1, SFTs1 );
130 REAL8 norm1 = sqrt( ( REAL8 )ret );
131
132 ret = scalarProductSFTVector( SFTs2, SFTs2 );
133 REAL8 norm2 = sqrt( ( REAL8 )ret );
134
135 ret = scalarProductSFTVector( SFTs1, SFTs2 );
136 REAL8 scalar = ( REAL8 ) ret;
137
138 ret = scalarProductSFTVector( diffs, diffs );
139 REAL8 normdiff = ( REAL8 ) ret;
140
141 REAL8 d1 = ( norm1 - norm2 ) / norm1;
142 maxd = fmax( maxd, d1 );
143 REAL8 d2 = 1.0 - scalar / ( norm1 * norm2 );
144 maxd = fmax( maxd, d2 );
145 REAL8 d3 = normdiff / ( norm1 * norm1 + norm2 * norm2 );
146 maxd = fmax( maxd, d3 );
147 REAL8 d4 = getMaxErrSFTVector( SFTs1, SFTs2 );
148 maxd = fmax( maxd, d4 );
149
150 if ( uvar.verbose ) {
151 printf( "TOTAL: (|x|-|y|)/|x|=%10.3e, 1-x.y/(|x||y|)=%10.3e, |x-y|^2/(|x|^2+|y|^2)=%10.3e, maxErrSFT=%10.3e\n", d1, d2, d3, d4 );
152 printf( "COMPARE: maxd=%10.3e, relErrMax=%10.3e\n", maxd, uvar.relErrorMax );
153 }
154
155 } /* combined total measures */
156
157
158 /* free memory */
159 XLALDestroySFTVector( SFTs1 );
160 XLALDestroySFTVector( SFTs2 );
161 XLALDestroySFTVector( diffs );
163
165
166 if ( maxd <= uvar.relErrorMax ) {
167 return 0;
168 } else {
169 if ( !uvar.quiet ) {
170 XLALPrintError( "Tolerance exceeded! maxd=%10.3e, relErrMax=%10.3e\n", maxd, uvar.relErrorMax );
171 }
172 return 1;
173 }
174
175} /* main */
176
177
178/*----------------------------------------------------------------------*/
179/* register all our "user-variables" */
180int
182{
183 XLAL_CHECK( uvar != NULL, XLAL_EINVAL );
184
185 /* set some defaults */
186 uvar->debug = lalDebugLevel;
187 uvar->verbose = 0;
188 uvar->quiet = 0;
189 uvar->relErrorMax = 1e-4;
190
191 /* now register all our user-variable */
192
193 XLAL_CHECK( XLALRegisterUvarMember( sftBname1, STRING, '1', REQUIRED, "Path and basefilename for SFTs1. Possibilities are:\n"
194 " - '<SFT file>;<SFT file>;...', where <SFT file> may contain wildcards\n - 'list:<file containing list of SFT files>'" ) == XLAL_SUCCESS, XLAL_EFUNC );
195 XLAL_CHECK( XLALRegisterUvarMember( sftBname2, STRING, '2', REQUIRED, "Path and basefilename for SFTs2. Possibilities are:\n"
196 " - '<SFT file>;<SFT file>;...', where <SFT file> may contain wildcards\n - 'list:<file containing list of SFT files>'" ) == XLAL_SUCCESS, XLAL_EFUNC );
197 XLAL_CHECK( XLALRegisterUvarMember( verbose, BOOLEAN, 'V', OPTIONAL, "Verbose output of differences" ) == XLAL_SUCCESS, XLAL_EFUNC );
198 XLAL_CHECK( XLALRegisterUvarMember( quiet, BOOLEAN, 'q', OPTIONAL, "No output of differences even on failure" ) == XLAL_SUCCESS, XLAL_EFUNC );
199 XLAL_CHECK( XLALRegisterUvarMember( relErrorMax, REAL8, 'e', OPTIONAL, "Maximal relative error acceptable to 'pass' comparison" ) == XLAL_SUCCESS, XLAL_EFUNC );
200
201 return XLAL_SUCCESS;
202
203} /* initUserVars() */
204
205/* for two SFTs: get maximal value of |X_k - Y_k|^2 / max(|X_k|^2,|Y_k|^2) */
206REAL4
207getMaxErrSFT( const SFTtype *sft1, const SFTtype *sft2 )
208{
209 UINT4 i;
210 REAL8 maxDiff, maxAmpl;
211 REAL4 maxErr;
212
213 maxDiff = 0;
214 maxAmpl = 0;
215 for ( i = 0; i < sft1->data->length; i++ ) {
216 REAL8 diff, A1, A2, Ampl;
217 REAL8 re1, re2, im1, im2;
218 re1 = crealf( sft1->data->data[i] );
219 im1 = cimagf( sft1->data->data[i] );
220 re2 = crealf( sft2->data->data[i] );
221 im2 = cimagf( sft2->data->data[i] );
222
223 diff = ( re1 - re2 ) * ( re1 - re2 ) + ( im1 - im2 ) * ( im1 - im2 );
224 A1 = re1 * re1 + im1 * im1;
225 A2 = re2 * re2 + im2 * im2;
226 Ampl = fmax( A1, A2 );
227
228 maxDiff = fmax( maxDiff, diff );
229 maxAmpl = fmax( maxAmpl, Ampl );
230
231 } /* for i */
232
233 maxErr = maxDiff / maxAmpl;
234
235 return ( maxErr );
236
237} /* getMaxErrSFT() */
238
239REAL4
240getMaxErrSFTVector( const SFTVector *sftvect1, const SFTVector *sftvect2 )
241{
242 UINT4 i;
243 REAL4 maxErr, thisErr;
244
245 maxErr = 0;
246
247 for ( i = 0; i < sftvect1->length; i++ ) {
248 thisErr = getMaxErrSFT( &( sftvect1->data[i] ), &( sftvect2->data[i] ) );
249 maxErr = fmax( maxErr, thisErr );
250 }
251
252 return maxErr;
253
254} /* getMaxErrSFTVector() */
255
256
257/*--------------------------------------------------
258 * implements a straightforward L2 scalar product of
259 * two time-series x_i and y_i : x*y = sum_i x_i y_i
260 * in Fourier-space, which is 2/N * Re( sum_i X_i Y*_i)
261 *--------------------------------------------------*/
262REAL4
263scalarProductSFT( const SFTtype *sft1, const SFTtype *sft2 )
264{
265 XLAL_CHECK_REAL4( ( sft1 != NULL ) && ( sft2 != NULL ), XLAL_EINVAL );
266 XLAL_CHECK_REAL4( sft1->data->length == sft2->data->length, XLAL_EINVAL );
267
268 /* we do the calculation in REAL8 to avoid accumulating roundoff-errors */
269 REAL8 prod = 0;
270 for ( UINT4 i = 0; i < sft1->data->length; i++ ) {
271 REAL8 xre = ( REAL8 )crealf( sft1->data->data[i] );
272 REAL8 xim = ( REAL8 )cimagf( sft1->data->data[i] );
273 REAL8 yre = ( REAL8 )crealf( sft2->data->data[i] );
274 REAL8 yim = ( REAL8 )cimagf( sft2->data->data[i] );
275
276 prod += xre * yre + xim * yim;
277
278 } /* for i < SFT-length */
279
280 prod *= 2.0 / ( ( REAL8 )sft1->data->length );
281
282 return ( REAL4 ) prod;
283
284} /* scalarProductSFT() */
285
286/*--------------------------------------------------
287 * extend the previous definition to an SFT-vector
288 * this is simply the sum of individual SFT-products
289 *--------------------------------------------------*/
290REAL4
291scalarProductSFTVector( const SFTVector *sftvect1, const SFTVector *sftvect2 )
292{
293 XLAL_CHECK_REAL4( ( sftvect1 != NULL ) && ( sftvect2 != NULL ), XLAL_EINVAL );
294 XLAL_CHECK_REAL4( sftvect1->length == sftvect2->length, XLAL_EINVAL );
295
296 REAL8 prod = 0;
297 for ( UINT4 i = 0; i < sftvect1->length; i++ ) {
298 REAL4 xy = scalarProductSFT( &( sftvect1->data[i] ), &( sftvect2->data[i] ) );
299 prod += ( REAL8 ) xy;
300 }
301
302 return ( REAL4 )prod;
303
304} /* scalarProductSFTVector() */
305
306
307/*--------------------------------------------------
308 * calculate the difference of two SFT-vectors
309 *--------------------------------------------------*/
310SFTVector *
311subtractSFTVectors( const SFTVector *sftvect1, const SFTVector *sftvect2 )
312{
313 XLAL_CHECK_NULL( ( sftvect1 != NULL ) && ( sftvect2 != NULL ), XLAL_EINVAL );
314 XLAL_CHECK_NULL( sftvect1->length == sftvect2->length, XLAL_EINVAL );
315
316 UINT4 numSFTs = sftvect1->length;
317 UINT4 numBins = sftvect1->data[0].data->length;
318
319 SFTVector *vect;
321
322 for ( UINT4 alpha = 0; alpha < numSFTs; alpha ++ ) {
323 for ( UINT4 j = 0; j < numBins; j++ ) {
324 vect->data[alpha].data->data[j] = sftvect1->data[alpha].data->data[j] - sftvect2->data[alpha].data->data[j];
325 }
326 } /* for alpha < numSFTs */
327
328 return vect;
329
330} /* subtractSFTVectors() */
int j
void LALCheckMemoryLeaks(void)
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
#define STRING(a)
double e
int main(int argc, char *argv[])
Definition: compareSFTs.c:56
int initUserVars(UserInput_t *uvar)
Definition: compareSFTs.c:181
REAL4 scalarProductSFTVector(const SFTVector *sftvect1, const SFTVector *sftvect2)
Definition: compareSFTs.c:291
REAL4 getMaxErrSFT(const SFTtype *sft1, const SFTtype *sft2)
Definition: compareSFTs.c:207
REAL4 scalarProductSFT(const SFTtype *sft1, const SFTtype *sft2)
Definition: compareSFTs.c:263
SFTVector * subtractSFTVectors(const SFTVector *sftvect1, const SFTVector *sftvect2)
Definition: compareSFTs.c:311
REAL4 getMaxErrSFTVector(const SFTVector *sftvect1, const SFTVector *sftvect2)
Definition: compareSFTs.c:240
char * XLALGPSToStr(char *s, const LIGOTimeGPS *t)
unsigned char BOOLEAN
double REAL8
#define XLAL_INIT_DECL(var,...)
char CHAR
uint32_t UINT4
int32_t INT4
float REAL4
#define lalDebugLevel
void XLALDestroySFTVector(SFTVector *vect)
XLAL interface to destroy an SFTVector.
Definition: SFTtypes.c:300
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
Definition: SFTcatalog.c:329
SFTVector * XLALLoadSFTs(const SFTCatalog *catalog, REAL8 fMin, REAL8 fMax)
Load the given frequency-band [fMin, fMax) (half-open) from the SFT-files listed in the SFT-'catalogu...
Definition: SFTfileIO.c:87
SFTCatalog * XLALSFTdataFind(const CHAR *file_pattern, const SFTConstraints *constraints)
Find the list of SFTs matching the file_pattern and satisfying the given constraints,...
Definition: SFTcatalog.c:71
SFTVector * XLALCreateSFTVector(UINT4 numSFTs, UINT4 numBins)
XLAL function to create an SFTVector of numSFT SFTs with SFTlen frequency-bins (which will be allocat...
Definition: SFTtypes.c:230
int XLALUserVarReadAllInput(BOOLEAN *should_exit, int argc, char *argv[], const LALVCSInfoList vcs_list)
void XLALDestroyUserVars(void)
#define XLALRegisterUvarMember(name, type, option, category,...)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_MAIN(assertion,...)
#define XLAL_CHECK_REAL4(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK_NULL(assertion,...)
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EINVAL
REAL8 XLALGPSDiff(const LIGOTimeGPS *t1, const LIGOTimeGPS *t0)
double alpha
CHAR name[LALNameLength]
COMPLEX8Sequence * data
A vector of COMPLEX8FrequencySeries.
COMPLEX8FrequencySeries * data
Pointer to the data array.
UINT4 length
Number of elements in array.
COMPLEX8 * data
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
Definition: SFTfileIO.h:238