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
ComputeFstatMCUpperLimit.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2008 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 * \author Karl Wette
22 * \file
23 * \ingroup lalpulsar_bin_Fstatistic
24 * \brief Computes an upper limit using Monte Carlo integration
25 * of the analytic F statistic signal model
26 */
27
28#include "config.h"
29
30#include <stdlib.h>
31#include <stdio.h>
32
33#include <gsl/gsl_math.h>
34#include <gsl/gsl_rng.h>
35#include <gsl/gsl_randist.h>
36#include <gsl/gsl_sf_bessel.h>
37#include <gsl/gsl_sf_gamma.h>
38#include <gsl/gsl_vector_int.h>
39#include <gsl/gsl_matrix.h>
40
41#include <lal/LALStdlib.h>
42#include <lal/LALString.h>
43#include <lal/LALError.h>
44#include <lal/LogPrintf.h>
45#include <lal/UserInput.h>
46#include <lal/SFTfileIO.h>
47#include <lal/LALInitBarycenter.h>
48#include <lal/DetectorStates.h>
49#include <lal/NormalizeSFTRngMed.h>
50#include <lal/ComputeFstat.h>
51#include <lal/ConfigFile.h>
52#include <lal/LALPulsarVCSInfo.h>
53
57REAL8 ran_ncx2_4( const gsl_rng *, REAL8 );
58gsl_vector_int *resize_histogram( gsl_vector_int *old_hist, size_t size );
59
60#define TRUE (1==1)
61#define FALSE (1==0)
62
63const REAL8 min_cosi = -1;
64const REAL8 max_cosi = 1;
67
68int main( int argc, char *argv[] )
69{
70
71 REAL8 alpha = 0.0;
72 REAL8 alpha_band = 0.0;
73 REAL8 delta = 0.0;
74 REAL8 delta_band = 0.0;
75 REAL8 freq = 0.0;
76 REAL8 freq_band = 0.0;
77 REAL8 twoFs = 0.0;
78 CHAR *mism_hist_file = NULL;
79 REAL8 max_mismatch = 0.0;
80 CHAR *sft_pattern = NULL;
81 CHAR *ephem_earth = NULL;
82 CHAR *ephem_sun = NULL;
84 REAL8 max_rel_err = 1.0e-3;
85 REAL8 h0_brake = 0.75;
86 INT4 MC_trials_init = 1e6;
87 REAL8 MC_trials_incr = 1.5;
88 INT4 MC_trials_reset_fac = 5e2;
89 REAL8 FDR = 0.05;
90 CHAR *output_file = NULL;
91 REAL8 twoF_pdf_hist_binw = 1.0;
92 CHAR *twoF_pdf_hist_file = NULL;
93 REAL8 initial_h0 = 0.0;
94
95 FILE *fp = NULL;
96 CHAR *cmdline = NULL;
97 gsl_matrix *mism_hist = NULL;
98 SFTCatalog *catalog = NULL;
99 MultiSFTVector *sfts = NULL;
100 EphemerisData *ephemeris = NULL;
101 MultiDetectorStateSeries *detector_states = NULL;
102 MultiPSDVector *rng_med = NULL;
103 MultiNoiseWeights *noise_weights = NULL;
104 BOOLEAN calc_ABC_coeffs = TRUE;
105 REAL8 A_coeff = 0.0;
106 REAL8 B_coeff = 0.0;
107 REAL8 C_coeff = 0.0;
108 gsl_rng *rng = NULL;
109 REAL8 h0 = 0.0;
110 REAL8 h0_prev = 0.0;
111 REAL8 dh0 = 0.0;
112 INT4 h0_iter = 0;
113 INT8 MC_trials = 0;
114 INT8 MC_iter = 0;
115 INT8 MC_trials_reset;
116 gsl_vector_int *twoF_pdf_hist = NULL;
117
118 ephem_earth = XLALStringDuplicate( "earth00-40-DE405.dat.gz" );
119 ephem_sun = XLALStringDuplicate( "sun00-40-DE405.dat.gz" );
120
121 /* Initialise LAL error handler */
123
124 /* Register command line arguments */
125 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &alpha, "alpha", REAL8, 'a', REQUIRED, "Right ascension in radians" ) == XLAL_SUCCESS, XLAL_EFUNC );
126 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &alpha, "alpha-band", REAL8, 'z', OPTIONAL, "Right ascension band" ) == XLAL_SUCCESS, XLAL_EFUNC );
127 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &delta, "delta", REAL8, 'd', REQUIRED, "Declination in radians" ) == XLAL_SUCCESS, XLAL_EFUNC );
128 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &delta, "delta-band", REAL8, 'c', OPTIONAL, "Declination band" ) == XLAL_SUCCESS, XLAL_EFUNC );
129 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &freq, "freq", REAL8, 'f', REQUIRED, "Starting frequency" ) == XLAL_SUCCESS, XLAL_EFUNC );
130 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &freq_band, "freq-band", REAL8, 'b', REQUIRED, "Frequency band" ) == XLAL_SUCCESS, XLAL_EFUNC );
131 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &twoFs, "loudest-2F", REAL8, 'F', REQUIRED, "Loudest 2F value in this band" ) == XLAL_SUCCESS, XLAL_EFUNC );
132 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &mism_hist_file, "mism-hist-file", STRING, 'M', OPTIONAL, "File containing the mismatch PDF histogram" ) == XLAL_SUCCESS, XLAL_EFUNC );
133 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &max_mismatch, "max-mismatch", REAL8, 'm', OPTIONAL, "Maximum mismatch to scale histogram to" ) == XLAL_SUCCESS, XLAL_EFUNC );
134 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &sft_pattern, "sft-patt", STRING, 'D', REQUIRED, "File pattern of the input SFTs. Possibilities are:\n"
135 " - '<SFT file>;<SFT file>;...', where <SFT file> may contain wildcards\n - 'list:<file containing list of SFT files>'" ) == XLAL_SUCCESS, XLAL_EFUNC );
136 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &ephem_earth, "ephem-earth", STRING, 'E', OPTIONAL, "Earth ephemeris file" ) == XLAL_SUCCESS, XLAL_EFUNC );
137 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &ephem_sun, "ephem-sun", STRING, 'S', OPTIONAL, "Sun ephemeris file" ) == XLAL_SUCCESS, XLAL_EFUNC );
138 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &rng_med_win, "rng-med-win", INT4, 'k', OPTIONAL, "Size of the running median window" ) == XLAL_SUCCESS, XLAL_EFUNC );
139 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &max_rel_err, "max-rel-err", REAL8, 'e', OPTIONAL, "Maximum error in h0 relative to previous value" ) == XLAL_SUCCESS, XLAL_EFUNC );
140 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &h0_brake, "h0-brake", REAL8, 'r', OPTIONAL, "h0 cannot change by more than this fraction of itself" ) == XLAL_SUCCESS, XLAL_EFUNC );
141 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &MC_trials_init, "init-MC-tri", INT4, 'I', OPTIONAL, "Initial number of MC int. trials" ) == XLAL_SUCCESS, XLAL_EFUNC );
142 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &MC_trials_incr, "MC-tri-incr", REAL8, 'i', OPTIONAL, "Multiply number of MC int. trials by this after each step" ) == XLAL_SUCCESS, XLAL_EFUNC );
143 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &MC_trials_reset_fac, "MC-tri-reset", INT4, 'Z', OPTIONAL, "Reset if no convergence after this ratio of initial MC int. trials" ) == XLAL_SUCCESS, XLAL_EFUNC );
144 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &FDR, "false-dism", REAL8, 'R', OPTIONAL, "Target false dismissal rate" ) == XLAL_SUCCESS, XLAL_EFUNC );
145 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &output_file, "output-file", STRING, 'o', OPTIONAL, "Output file for the upper limit and other info (defaults to stdout)" ) == XLAL_SUCCESS, XLAL_EFUNC );
146 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &twoF_pdf_hist_binw, "2F-pdf-hist-binw", REAL8, 'B', OPTIONAL, "Bin width of the histogram of the non-central 2F distribution" ) == XLAL_SUCCESS, XLAL_EFUNC );
147 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &twoF_pdf_hist_file, "2F-pdf-hist-file", STRING, 'H', OPTIONAL, "Output file for the histogram of the non-central 2F distribution" ) == XLAL_SUCCESS, XLAL_EFUNC );
148 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &initial_h0, "initial-h0", REAL8, 0, DEVELOPER, "Initial guess of h0 (default: automatic)" ) == XLAL_SUCCESS, XLAL_EFUNC );
149
150 /* Get command line arguments */
151 BOOLEAN should_exit = 0;
153 if ( should_exit ) {
154 return EXIT_FAILURE;
155 }
156
157 /* Load the mismatch PDF histogram if available */
158 if ( XLALUserVarWasSet( &mism_hist_file ) ) {
159
160 size_t i;
161 LALParsedDataFile *file = NULL;
162
163 /* Check the mismatch */
164 if ( max_mismatch <= 0.0 ) {
165 XLALPrintError( "Maximum mismatch must be strictly positive\n" );
166 return EXIT_FAILURE;
167 }
168
169 /* Load the file */
170 XLAL_CHECK( XLALParseDataFile( &file, mism_hist_file ) == XLAL_SUCCESS, XLAL_EINVAL );
171
172 /* Allocate memory */
173 if ( ( mism_hist = gsl_matrix_alloc( file->lines->nTokens, 3 ) ) == NULL ) {
174 XLALPrintError( "Couldn't allocate a gsl_matrix\n" );
175 return EXIT_FAILURE;
176 }
177
178 /* Read in histogram */
179 for ( i = 0; i < mism_hist->size1; ++i ) {
180 double left_bin, right_bin, prob_bin;
181 if ( sscanf( file->lines->tokens[i], "%lf %lf %lf", &left_bin, &right_bin, &prob_bin ) != 3 ) {
182 XLALPrintError( "Couldn't parse line %zu of '%s'\n", i, mism_hist_file );
183 return EXIT_FAILURE;
184 }
185 if ( left_bin >= right_bin || prob_bin < 0.0 ) {
186 XLALPrintError( "Invalid syntax: line %zu of '%s'\n", i, mism_hist_file );
187 return EXIT_FAILURE;
188 }
189 gsl_matrix_set( mism_hist, i, 0, left_bin );
190 gsl_matrix_set( mism_hist, i, 1, right_bin );
191 gsl_matrix_set( mism_hist, i, 2, prob_bin );
192 }
193
194 /* Cleanup */
196
197 /* Rescale histogram to maximum mismatch */
198 {
199 double max_bin = 0.0;
200 for ( i = 0; i < mism_hist->size1; ++i ) {
201 const double right_bin = gsl_matrix_get( mism_hist, i, 1 );
202 if ( max_bin < right_bin ) {
203 max_bin = right_bin;
204 }
205 }
206 for ( i = 0; i < mism_hist->size1; ++i ) {
207 gsl_matrix_set( mism_hist, i, 0, gsl_matrix_get( mism_hist, i, 0 ) / max_bin * max_mismatch );
208 gsl_matrix_set( mism_hist, i, 1, gsl_matrix_get( mism_hist, i, 1 ) / max_bin * max_mismatch );
209 }
210 }
211
212 /* Normalise histogram to unit area */
213 {
214 double total_area = 0.0;
215 for ( i = 0; i < mism_hist->size1; ++i ) {
216 const double left_bin = gsl_matrix_get( mism_hist, i, 0 );
217 const double right_bin = gsl_matrix_get( mism_hist, i, 1 );
218 const double prob_bin = gsl_matrix_get( mism_hist, i, 2 );
219 total_area += prob_bin * ( right_bin - left_bin );
220 }
221 for ( i = 0; i < mism_hist->size1; ++i ) {
222 gsl_matrix_set( mism_hist, i, 2, gsl_matrix_get( mism_hist, i, 2 ) / total_area );
223 }
224 }
225
226 }
227
228 /* Load the SFTs */
229 {
230 SFTConstraints XLAL_INIT_DECL( constraints );
231 REAL8 extra = 0.0, f_min = 0.0, f_max = 0.0;
232
233 /* Load the catalog */
234 LogPrintf( LOG_DEBUG, "Loading SFT catalog ... " );
235 XLAL_CHECK_MAIN( ( catalog = XLALSFTdataFind( sft_pattern, &constraints ) ) != NULL, XLAL_EFUNC );
236 if ( !catalog || catalog->length == 0 ) {
237 XLALPrintError( "Couldn't find SFTs matching '%s'\n", sft_pattern );
238 return EXIT_FAILURE;
239 }
240 LogPrintfVerbatim( LOG_DEBUG, "done: %i SFTs starting at GPS %i\n",
241 catalog->length, catalog->data[0].header.epoch.gpsSeconds );
242
243 /* Determine the frequency range */
244 extra = catalog->data[0].header.deltaF * ( rng_med_win / 2 + 1 );
245 f_min = freq - extra;
246 f_max = freq + extra + freq_band;
247
248 /* Load the SFTs */
249 LogPrintf( LOG_DEBUG, "Loading SFTs (%f to %f) ... ", f_min, f_max );
250 XLAL_CHECK_MAIN( ( sfts = XLALLoadMultiSFTs( catalog, f_min, f_max ) ) != NULL, XLAL_EFUNC );
251 LogPrintfVerbatim( LOG_DEBUG, "done\n" );
252 }
253
254 /* Load the ephemeris data */
255 LogPrintf( LOG_DEBUG, "Loading ephemeris ... " );
256 XLAL_CHECK_MAIN( ( ephemeris = XLALInitBarycenter( ephem_earth, ephem_sun ) ) != NULL, XLAL_EFUNC );
257 LogPrintfVerbatim( LOG_DEBUG, "done\n" );
258
259 /* Get the detector states */
260 LogPrintf( LOG_DEBUG, "Calculating detector states ... " );
261 const REAL8 tOffset = 0.5 / sfts->data[0]->data[0].deltaF;
262 XLAL_CHECK_MAIN( ( detector_states = XLALGetMultiDetectorStatesFromMultiSFTs( sfts, ephemeris, tOffset ) ) != NULL, XLAL_EFUNC );
263 LogPrintfVerbatim( LOG_DEBUG, "done\n" );
264
265 /* Normalise SFTs and compute noise weights */
266 LogPrintf( LOG_DEBUG, "Normalising SFTs and computing noise weights ... " );
267 XLAL_CHECK_MAIN( ( rng_med = XLALNormalizeMultiSFTVect( sfts, rng_med_win, NULL ) ) != NULL, XLAL_EFUNC );
268 XLAL_CHECK_MAIN( ( noise_weights = XLALComputeMultiNoiseWeights( rng_med, rng_med_win, 0 ) ) != NULL, XLAL_EFUNC );
269 LogPrintfVerbatim( LOG_DEBUG, "done\n" );
270
271 /* Cleanup */
272 XLALDestroySFTCatalog( catalog );
274 XLALDestroyEphemerisData( ephemeris );
275 XLALDestroyMultiPSDVector( rng_med );
276
277 /* Initialise the random number generator */
278 {
279 unsigned long seed;
280 FILE *fpr = NULL;
281
282 if ( ( rng = gsl_rng_alloc( gsl_rng_mt19937 ) ) == NULL ) {
283 XLALPrintError( "Couldn't allocate a gsl_rng\n" );
284 return EXIT_FAILURE;
285 }
286 /*
287 * Note: /dev/random can be slow after the first few accesses, which is why we're using urandom instead.
288 * [Cryptographic safety isn't a concern here at all]
289 */
290 if ( ( fpr = fopen( "/dev/urandom", "r" ) ) == NULL ) {
291 XLALPrintError( "Couldn't open '/dev/urandom'\n" );
292 return EXIT_FAILURE;
293 }
294 if ( fread( &seed, sizeof( seed ), 1, fpr ) != 1 ) {
295 XLALPrintError( "Couldn't read from '/dev/urandom'\n" );
296 return EXIT_FAILURE;
297 }
298 fclose( fpr );
299 gsl_rng_set( rng, seed );
300 }
301
302 /* Calculate the AM coefficients at least once */
303 calc_ABC_coeffs = calc_AM_coeffs( rng,
304 alpha, alpha_band,
305 delta, delta_band,
306 detector_states, noise_weights,
307 &A_coeff, &B_coeff, &C_coeff );
308
309 /* Begin iterations to find h0 */
310 MC_trials_reset = ( ( INT8 )MC_trials_reset_fac * ( ( INT8 )MC_trials_init ) );
311 do {
312
313 /* Open the output file */
314 if ( XLALUserVarWasSet( &output_file ) ) {
315 if ( ( fp = fopen( output_file, "wb" ) ) == NULL ) {
316 XLALPrintError( "Couldn't open output file '%s'\n", output_file );
317 return EXIT_FAILURE;
318 }
319 } else {
320 fp = stdout;
321 }
323 /** \deprecated FIXME: the following code uses obsolete CVS ID tags.
324 * It should be modified to use git version information. */
325 fprintf( fp, "%%%% %s\n%%%% %s\n", "$Id$", cmdline );
326 LALFree( cmdline );
327 cmdline = NULL;
328 fprintf( fp, "alpha=%0.4f alpha_band=%0.4f\n", alpha, alpha_band );
329 fprintf( fp, "delta=%0.4f delta_band=%0.4f\n", delta, delta_band );
330 fprintf( fp, "freq=%0.4f freq_band=%0.4f\n", freq, freq_band );
331
332 /* Use initial h0 guess if given (first time only) */
333 if ( initial_h0 > 0.0 ) {
334 h0 = initial_h0;
335 initial_h0 = 0.0;
336 }
337 /* Otherwise compute first guess at h0 */
338 else {
339 h0 = twoFs * pow( A_coeff * B_coeff - C_coeff * C_coeff, -0.25 );
340 h0 *= GSL_MAX( 0.5, 1.0 + gsl_ran_gaussian( rng, 0.1 ) );
341 }
342
343 /* Begin Newton-Raphson iteration to find h0 */
344 h0_prev = GSL_POSINF;
345#define H0_ERROR fabs(1.0 - (h0 / h0_prev))
346 h0_iter = 0;
347 MC_trials = MC_trials_init;
348 while ( TRUE ) {
349
350 /* Integrand and its derivative w.r.t. h0 */
351 REAL8 J = 0.0;
352 REAL8 dJ = 0.0;
353
354 INT4 twoF_pdf_FD = 0;
355 REAL8 twoF_pdf_FDR = 0.0;
356
357 /* Output at beginning of loop */
358 LogPrintf( LOG_DEBUG, "Beginning h0 loop %2i with h0=%0.5e, dh0=% 0.5e, MC_trials=%" LAL_INT8_FORMAT "\n", h0_iter, h0, dh0, MC_trials );
359 fprintf( fp, "MC_trials=%" LAL_INT8_FORMAT, MC_trials );
360
361 /* Destroy any previous histogram */
362 if ( twoF_pdf_hist ) {
363 gsl_vector_int_free( twoF_pdf_hist );
364 twoF_pdf_hist = NULL;
365 }
366
367 /* Begin Monte Carlo integration to find J and dJ */
368 MC_iter = MC_trials;
369 while ( MC_iter-- ) {
370
371 size_t i;
372
373 REAL8 cosi = 0.0;
374 REAL8 psi = 0.0;
375 REAL8 twoF = 0.0;
376
377 REAL8 A1 = 0.0;
378 REAL8 A2 = 0.0;
379 REAL8 A3 = 0.0;
380 REAL8 A4 = 0.0;
381
382 REAL8 rho2 = 0.0;
383 REAL8 mism_rho2 = 0.0;
384
385 REAL8 mismatch = 0.0;
386 REAL8 prob_mismatch = 1.0;
387
388 REAL8 mismatch_from_pdf = 0.0;
389 REAL8 twoF_from_pdf = 0.0;
390
391 /* Generate random cosi, psi, and twoF */
392 cosi = gsl_ran_flat( rng, min_cosi, max_cosi );
393 psi = gsl_ran_flat( rng, min_psi, max_psi );
394 twoF = gsl_ran_flat( rng, 0, twoFs );
395
396 /* Calculate the AM coefficients if needed */
397 if ( calc_ABC_coeffs )
398 calc_ABC_coeffs = calc_AM_coeffs( rng,
399 alpha, alpha_band,
400 delta, delta_band,
401 detector_states, noise_weights,
402 &A_coeff, &B_coeff, &C_coeff );
403
404 /* Compute the amplitude coefficients vector A */
405 A1 = 0.5 * h0 * ( 1 + cosi * cosi ) * cos( 2 * psi );
406 A2 = 0.5 * h0 * ( 1 + cosi * cosi ) * sin( 2 * psi );
407 A3 = -1.0 * h0 * cosi * sin( 2 * psi );
408 A4 = 1.0 * h0 * cosi * cos( 2 * psi );
409
410 /* Compute the optimal signal to noise ratio rho^2 */
411 rho2 = ( A1 * A1 * A_coeff + A2 * A2 * B_coeff +
412 A1 * A2 * C_coeff + A2 * A1 * C_coeff +
413 A3 * A3 * A_coeff + A4 * A4 * B_coeff +
414 A3 * A4 * C_coeff + A4 * A3 * C_coeff );
415
416 /* Generate random mismatch */
417 if ( XLALUserVarWasSet( &mism_hist_file ) ) {
418
419 mismatch = gsl_ran_flat( rng, 0.0, max_mismatch );
420
421 prob_mismatch = 0.0;
422 for ( i = 0; i < mism_hist->size1; ++i ) {
423 const double left_bin = gsl_matrix_get( mism_hist, i, 0 );
424 const double right_bin = gsl_matrix_get( mism_hist, i, 1 );
425 const double prob_bin = gsl_matrix_get( mism_hist, i, 2 );
426 if ( left_bin <= mismatch && mismatch < right_bin ) {
427 prob_mismatch = prob_bin;
428 break;
429 }
430 }
431
432 }
433 mism_rho2 = ( 1.0 - mismatch ) * rho2;
434
435 /* Add to the integrand and its derivative w.r.t. h0 */
436 J += pdf_ncx2_4( mism_rho2, twoF ) * prob_mismatch;
437 dJ += 2.0 * mism_rho2 / h0 * d_pdf_ncx2_4( mism_rho2, twoF ) * prob_mismatch;
438
439 /* Break if J and dJ failed */
440 if ( gsl_isnan( J ) || gsl_isnan( dJ ) ) {
441 break;
442 }
443
444 /* Draw a random mismatch from the mismatch histogram */
445 if ( XLALUserVarWasSet( &mism_hist_file ) ) {
446
447 REAL8 mism_pdf_cumul_prob_bin = gsl_ran_flat( rng, 0.0, 1.0 );
448
449 mismatch_from_pdf = max_mismatch;
450 for ( i = 0; i < mism_hist->size1; ++i ) {
451 const double left_bin = gsl_matrix_get( mism_hist, i, 0 );
452 const double right_bin = gsl_matrix_get( mism_hist, i, 1 );
453 const double prob_bin = gsl_matrix_get( mism_hist, i, 2 );
454 if ( mism_pdf_cumul_prob_bin < ( right_bin - left_bin ) ) {
455 mismatch_from_pdf = left_bin + mism_pdf_cumul_prob_bin / prob_bin;
456 break;
457 }
458 mism_pdf_cumul_prob_bin -= prob_bin * ( right_bin - left_bin );
459 }
460
461 }
462 mism_rho2 = ( 1.0 - mismatch_from_pdf ) * rho2;
463
464 /* Draw a 2F value from the non-central chi-square with parameter rho^2 */
465 twoF_from_pdf = ran_ncx2_4( rng, mism_rho2 );
466
467 /* Count 2F value if it is below threshold */
468 if ( twoF_from_pdf < twoFs ) {
469 ++twoF_pdf_FD;
470 }
471
472 /* Add 2F to histogram if needed */
473 if ( XLALUserVarWasSet( &twoF_pdf_hist_file ) ) {
474
475 /* Compute bin */
476 const size_t bin = twoF_from_pdf / twoF_pdf_hist_binw;
477
478 /* Resize histogram vector if needed */
479 if ( !twoF_pdf_hist || bin >= twoF_pdf_hist->size )
480 if ( NULL == ( twoF_pdf_hist = resize_histogram( twoF_pdf_hist, bin + 1 ) ) ) {
481 XLALPrintError( "\nCouldn't (re)allocate 'twoF_pdf_hist'\n" );
482 return EXIT_FAILURE;
483 }
484
485 /* Add to bin */
486 gsl_vector_int_set( twoF_pdf_hist, bin,
487 gsl_vector_int_get( twoF_pdf_hist, bin ) + 1 );
488
489 }
490
491 }
492
493 /* If J and dJ failed, reduce h0 and try again */
494 if ( gsl_isnan( J ) || gsl_isnan( dJ ) ) {
495 h0 /= 2.0;
496 LogPrintf( LOG_DEBUG, "Reducing h0 to %0.4e and starting again because J=%0.4e, dJ=%0.4e\n", h0, J, dJ );
497 h0_iter = 0;
498 MC_trials = MC_trials_init;
499 continue;
500 }
501
502 /* Monte Carlo normalisation: integrate over 2F,mismatch but average cosi,psi; divide by number of trials */
503 {
504 REAL8 MC_norm = twoFs / MC_trials;
505 if ( XLALUserVarWasSet( &mism_hist_file ) ) {
506 MC_norm *= max_mismatch;
507 }
508
509 J *= MC_norm;
510 dJ *= MC_norm;
511 }
512
513 /* Compute the false dismissal rate from 2F distribution */
514 twoF_pdf_FDR = ( 1.0 * twoF_pdf_FD ) / MC_trials;
515
516 /* Compute the increment in h0 from Newton-Raphson */
517 dh0 = ( FDR - J ) / dJ;
518
519 /* Limit the increment in h0 to |h0 * h0_brake| */
520 dh0 = GSL_SIGN( dh0 ) * GSL_MIN( fabs( dh0 ), fabs( h0 * h0_brake ) );
521
522 /* Output at end of loop */
523 fprintf( fp, " h0=%0.5e FDR_MC_int=%0.4f FDR_2F_dist=%0.4f\n", h0, J, twoF_pdf_FDR );
524 fflush( fp );
525 LogPrintf( LOG_DEBUG, "Ending h0 loop %2i with error=%0.4e, FDR(MC int.)=%0.4f, FDR(2F dist.)=%0.4f\n", h0_iter, H0_ERROR, J, twoF_pdf_FDR );
526
527 /* Done if error condition or maximum number of trials */
528 if ( H0_ERROR <= max_rel_err || MC_trials >= MC_trials_reset ) {
529 break;
530 }
531
532 /* Increment h0 */
533 h0_prev = h0;
534 h0 += dh0;
535
536 /* Increase iteration count and number of MC trials */
537 ++h0_iter;
538 MC_trials = ( INT8 )ceil( MC_trials * MC_trials_incr );
539
540 }
541#undef H0_ERROR
542
543 /* Close the output file */
544 if ( XLALUserVarWasSet( &output_file ) ) {
545 fprintf( fp, "%%DONE\n" );
546 fclose( fp );
547 }
548
549 /* If number of MC trials exceeded reset */
550 if ( MC_trials >= MC_trials_reset ) {
551 LogPrintf( LOG_DEBUG, "Failed to converge after %i iterations (MC_trials=%" LAL_INT8_FORMAT "): trying again ...\n", h0_iter, MC_trials );
552 }
553
554 } while ( MC_trials >= MC_trials_reset );
555
556 /* Write 2F histogram if needed */
557 if ( XLALUserVarWasSet( &twoF_pdf_hist_file ) ) {
558
559 size_t i;
560 FILE *fpH;
561
562 if ( ( fpH = fopen( twoF_pdf_hist_file, "wb" ) ) == NULL ) {
563 XLALPrintError( "Couldn't open histogram file '%s'\n", twoF_pdf_hist_file );
564 return EXIT_FAILURE;
565 }
567 /** \deprecated FIXME: the following code uses obsolete CVS ID tags.
568 * It should be modified to use git version information. */
569 fprintf( fpH, "%%%% %s\n%%%% %s\n", "$Id$", cmdline );
570 LALFree( cmdline );
571 cmdline = NULL;
572
573 for ( i = 0; i < twoF_pdf_hist->size; ++i )
574 fprintf( fpH, "%0.3g %0.3g %i\n",
575 twoF_pdf_hist_binw * i,
576 twoF_pdf_hist_binw * ( i + 1 ),
577 gsl_vector_int_get( twoF_pdf_hist, i ) );
578
579 fprintf( fpH, "%%DONE\n" );
580 fclose( fpH );
581
582 }
583
584 /* Cleanup */
586 if ( mism_hist ) {
587 gsl_matrix_free( mism_hist );
588 }
589 XLALDestroyMultiDetectorStateSeries( detector_states );
590 XLALDestroyMultiNoiseWeights( noise_weights );
591 gsl_rng_free( rng );
592 if ( twoF_pdf_hist ) {
593 gsl_vector_int_free( twoF_pdf_hist );
594 }
596
597 return EXIT_SUCCESS;
598
599}
600
601/* Compute the AM coefficients */
603 gsl_rng *rng,
604 REAL8 alpha,
605 REAL8 alpha_band,
606 REAL8 delta,
607 REAL8 delta_band,
608 MultiDetectorStateSeries *detector_states,
609 MultiNoiseWeights *noise_weights,
610 REAL8 *A_coeff,
611 REAL8 *B_coeff,
612 REAL8 *C_coeff )
613{
614
615 MultiAMCoeffs *AM_coeffs = NULL;
617
618 /* Generate a random sky position */
619 sky.system = COORDINATESYSTEM_EQUATORIAL;
620 sky.longitude = gsl_ran_flat( rng, alpha, alpha + alpha_band );
621 sky.latitude = gsl_ran_flat( rng, delta, delta + delta_band );
622
623 /* Calculate and noise-weigh the AM coefficients */
624 XLAL_CHECK_MAIN( ( AM_coeffs = XLALComputeMultiAMCoeffs( detector_states, noise_weights, sky ) ) != NULL, XLAL_EFUNC );
625 *A_coeff = AM_coeffs->Mmunu.Ad * AM_coeffs->Mmunu.Sinv_Tsft;
626 *B_coeff = AM_coeffs->Mmunu.Bd * AM_coeffs->Mmunu.Sinv_Tsft;
627 *C_coeff = AM_coeffs->Mmunu.Cd * AM_coeffs->Mmunu.Sinv_Tsft;
628
629 /* Correct for use of DOUBLE-SIDED PSD by AM coefficient functions */
630 /* *A_coeff *= 0.5; *B_coeff *= 0.5; *C_coeff *= 0.5;
631 RP: commented-out as the field Sinv_Tsft has now refers to the single-sided PSD!
632 */
633
634 /* Cleanup */
635 XLALDestroyMultiAMCoeffs( AM_coeffs );
636
637 /* Return if AM coefficients need to be calculated again */
638 return ( alpha_band != 0.0 || delta_band != 0.0 );
639
640}
641
642/* PDF of non-central chi-square with 4 degrees
643 of freedom and non-centrality parameter lambda */
645{
646
647 const REAL8 z = sqrt( x * lambda );
648
649 gsl_sf_result I1;
650 gsl_error_handler_t *h = NULL;
651
652 /* Compute the Bessel functions */
653 h = gsl_set_error_handler_off();
654 if ( gsl_sf_bessel_In_e( 1, z, &I1 ) != GSL_SUCCESS ) {
655 gsl_set_error_handler( h );
656 return GSL_NAN;
657 }
658 gsl_set_error_handler( h );
659
660 /* Compute the PDF */
661 return 0.5 * exp( -0.5 * ( x + lambda ) ) * sqrt( x / lambda ) * I1.val;
662
663}
664
665/* Derivative of the PDF of non-central chi-square with 4 degrees
666 of freedom and non-centrality parameter lambda w.r.t. lambda */
668{
669
670 const REAL8 z = sqrt( x * lambda );
671
672 int i;
673 gsl_sf_result In[3];
674 gsl_error_handler_t *h = NULL;
675
676 /* Compute the Bessel functions */
677 h = gsl_set_error_handler_off();
678 for ( i = 0; i < 3; ++i )
679 if ( gsl_sf_bessel_In_e( i, z, &In[i] ) != GSL_SUCCESS ) {
680 gsl_set_error_handler( h );
681 return GSL_NAN;
682 }
683 gsl_set_error_handler( h );
684
685 /* Compute the derivative of the PDF */
686 return 0.25 * exp( -0.5 * ( x + lambda ) ) * ( 0.5 * x / lambda * ( In[0].val + In[2].val ) - sqrt( x / lambda ) * ( 1 + 1 / lambda ) * In[1].val );
687
688}
689
690/* Random number drawn from a non-central chi-square with 4 degrees
691 of freedom and non-centrality parameter lambda */
692REAL8 ran_ncx2_4( const gsl_rng *rng, REAL8 lambda )
693{
694
695 const REAL8 a = sqrt( lambda / 4 );
696
697 int i;
698 REAL8 x = 0.0;
699
700 for ( i = 0; i < 4; ++i ) {
701 x += pow( gsl_ran_gaussian( rng, 1.0 ) + a, 2.0 );
702 }
703
704 return x;
705
706}
707
708/* Resize histogram */
709gsl_vector_int *resize_histogram( gsl_vector_int *old_hist, size_t size )
710{
711 gsl_vector_int *new_hist = gsl_vector_int_alloc( size );
712 XLAL_CHECK_NULL( new_hist != NULL, XLAL_ENOMEM );
713 gsl_vector_int_set_zero( new_hist );
714 if ( old_hist != NULL ) {
715 gsl_vector_int_view old = gsl_vector_int_subvector( old_hist, 0, GSL_MIN( old_hist->size, size ) );
716 gsl_vector_int_view new = gsl_vector_int_subvector( new_hist, 0, GSL_MIN( old_hist->size, size ) );
717 gsl_vector_int_memcpy( &new.vector, &old.vector );
718 gsl_vector_int_free( old_hist );
719 }
720 return new_hist;
721}
int main(int argc, char *argv[])
REAL8 ran_ncx2_4(const gsl_rng *, REAL8)
#define H0_ERROR
const REAL8 min_psi
const REAL8 min_cosi
const REAL8 max_cosi
BOOLEAN calc_AM_coeffs(gsl_rng *, REAL8, REAL8, REAL8, REAL8, MultiDetectorStateSeries *, MultiNoiseWeights *, REAL8 *, REAL8 *, REAL8 *)
REAL8 pdf_ncx2_4(REAL8, REAL8)
gsl_vector_int * resize_histogram(gsl_vector_int *old_hist, size_t size)
#define TRUE
const REAL8 max_psi
REAL8 d_pdf_ncx2_4(REAL8, REAL8)
lal_errhandler_t lal_errhandler
int LAL_ERR_EXIT(LALStatus *stat, const char *func, const char *file, const int line, volatile const char *id)
void LALCheckMemoryLeaks(void)
#define LALFree(p)
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
static double double delta
#define STRING(a)
#define fprintf
const double rho2
const FstatOptionalArgs FstatOptionalArgsDefaults
Global initializer for setting FstatOptionalArgs to default values.
Definition: ComputeFstat.c:90
int XLALParseDataFile(LALParsedDataFile **cfgdata, const CHAR *fname)
void XLALDestroyParsedDataFile(LALParsedDataFile *cfgdata)
void XLALDestroyMultiDetectorStateSeries(MultiDetectorStateSeries *mdetStates)
Helper function to get rid of a multi-IFO DetectorStateSeries Note, this is "NULL-robust" in the sens...
MultiDetectorStateSeries * XLALGetMultiDetectorStatesFromMultiSFTs(const MultiSFTVector *multiSFTs, const EphemerisData *edat, REAL8 tOffset)
Get the 'detector state' (ie detector-tensor, position, velocity, etc) for the given multi-vector of ...
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
void XLALDestroyMultiAMCoeffs(MultiAMCoeffs *multiAMcoef)
Destroy a MultiAMCoeffs structure.
Definition: LALComputeAM.c:469
MultiAMCoeffs * XLALComputeMultiAMCoeffs(const MultiDetectorStateSeries *multiDetStates, const MultiNoiseWeights *multiWeights, SkyPosition skypos)
Multi-IFO version of XLALComputeAMCoeffs().
Definition: LALComputeAM.c:379
#define LAL_PI_4
unsigned char BOOLEAN
double REAL8
#define XLAL_INIT_DECL(var,...)
int64_t INT8
char CHAR
int32_t INT4
#define LAL_INT8_FORMAT
char char * XLALStringDuplicate(const char *s)
void void LogPrintfVerbatim(LogLevel_t, const char *format,...) _LAL_GCC_PRINTF_FORMAT_(2
void LogPrintf(LogLevel_t, const char *format,...) _LAL_GCC_PRINTF_FORMAT_(2
LOG_DEBUG
MultiPSDVector * XLALNormalizeMultiSFTVect(MultiSFTVector *multsft, UINT4 blockSize, const MultiNoiseFloor *assumeSqrtSX)
Function for normalizing a multi vector of SFTs in a multi IFO search and returns the running-median ...
void XLALDestroyMultiPSDVector(MultiPSDVector *multvect)
Destroy a multi PSD-vector.
Definition: PSDutils.c:102
MultiNoiseWeights * XLALComputeMultiNoiseWeights(const MultiPSDVector *rngmed, UINT4 blocksRngMed, UINT4 excludePercentile)
Computes weight factors arising from MultiSFTs with different noise floors.
Definition: PSDutils.c:285
void XLALDestroyMultiNoiseWeights(MultiNoiseWeights *weights)
Destroy a MultiNoiseWeights object.
Definition: PSDutils.c:172
static const INT4 a
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
Definition: SFTcatalog.c:329
MultiSFTVector * XLALLoadMultiSFTs(const SFTCatalog *inputCatalog, REAL8 fMin, REAL8 fMax)
Function to load a catalog of SFTs from possibly different detectors.
Definition: SFTfileIO.c:416
void XLALDestroyMultiSFTVector(MultiSFTVector *multvect)
Destroy a multi SFT-vector.
Definition: SFTtypes.c:424
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
COORDINATESYSTEM_EQUATORIAL
int XLALUserVarReadAllInput(BOOLEAN *should_exit, int argc, char *argv[], const LALVCSInfoList vcs_list)
void XLALDestroyUserVars(void)
void CHAR * XLALUserVarGetLog(UserVarLogFormat format)
#define XLALRegisterNamedUvar(cvar, name, type, option, category,...)
int XLALUserVarWasSet(const void *cvar)
UVAR_LOGFMT_CMDLINE
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_MAIN(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK_NULL(assertion,...)
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EINVAL
size
output_file
double alpha
REAL8 Sinv_Tsft
normalization-factor (using single-sided PSD!)
Definition: LALComputeAM.h:133
COMPLEX8FrequencySeries * data
Pointer to the data array.
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
UINT4 runningMedianWindow
If SFT noise weights are calculated from the SFTs, the running median window length to use.
Definition: ComputeFstat.h:141
Multi-IFO container for antenna-pattern coefficients and atenna-pattern matrix .
Definition: LALComputeAM.h:137
AntennaPatternMatrix Mmunu
antenna-pattern matrix
Definition: LALComputeAM.h:143
Multi-IFO time-series of DetectorStates.
One noise-weight (number) per SFT (therefore indexed over IFOs and SFTs.
Definition: PSDutils.h:71
A collection of PSD vectors – one for each IFO in a multi-IFO search.
Definition: PSDutils.h:62
A collection of SFT vectors – one for each IFO in a multi-IFO search.
Definition: SFTfileIO.h:179
SFTVector ** data
sftvector for each ifo
Definition: SFTfileIO.h:184
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
Definition: SFTfileIO.h:238
SFTDescriptor * data
array of data-entries describing matched SFTs
Definition: SFTfileIO.h:243
UINT4 length
number of SFTs in catalog
Definition: SFTfileIO.h:242
'Constraints' for SFT-matching: which detector, within which time-stretch and which timestamps exactl...
Definition: SFTfileIO.h:212
SFTtype header
SFT-header info.
Definition: SFTfileIO.h:228
double f_min
double f_max