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>
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>
68int main(
int argc,
char *argv[] )
72 REAL8 alpha_band = 0.0;
74 REAL8 delta_band = 0.0;
76 REAL8 freq_band = 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;
91 REAL8 twoF_pdf_hist_binw = 1.0;
92 CHAR *twoF_pdf_hist_file = NULL;
93 REAL8 initial_h0 = 0.0;
97 gsl_matrix *mism_hist = NULL;
115 INT8 MC_trials_reset;
116 gsl_vector_int *twoF_pdf_hist = NULL;
135 " - '<SFT file>;<SFT file>;...', where <SFT file> may contain wildcards\n - 'list:<file containing list of SFT files>'" ) ==
XLAL_SUCCESS,
XLAL_EFUNC );
164 if ( max_mismatch <= 0.0 ) {
173 if ( ( mism_hist = gsl_matrix_alloc(
file->lines->nTokens, 3 ) ) == NULL ) {
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 );
185 if ( left_bin >= right_bin || prob_bin < 0.0 ) {
186 XLALPrintError(
"Invalid syntax: line %zu of '%s'\n",
i, mism_hist_file );
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 );
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 ) {
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 );
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 );
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 );
236 if ( !catalog || catalog->
length == 0 ) {
237 XLALPrintError(
"Couldn't find SFTs matching '%s'\n", sft_pattern );
282 if ( ( rng = gsl_rng_alloc( gsl_rng_mt19937 ) ) == NULL ) {
290 if ( ( fpr = fopen(
"/dev/urandom",
"r" ) ) == NULL ) {
294 if ( fread( &
seed,
sizeof(
seed ), 1, fpr ) != 1 ) {
299 gsl_rng_set( rng,
seed );
306 detector_states, noise_weights,
307 &A_coeff, &B_coeff, &C_coeff );
310 MC_trials_reset = ( (
INT8 )MC_trials_reset_fac * ( (
INT8 )MC_trials_init ) );
325 fprintf(
fp,
"%%%% %s\n%%%% %s\n",
"$Id$", cmdline );
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 );
333 if ( initial_h0 > 0.0 ) {
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 ) );
344 h0_prev = GSL_POSINF;
345#define H0_ERROR fabs(1.0 - (h0 / h0_prev))
347 MC_trials = MC_trials_init;
354 INT4 twoF_pdf_FD = 0;
355 REAL8 twoF_pdf_FDR = 0.0;
362 if ( twoF_pdf_hist ) {
363 gsl_vector_int_free( twoF_pdf_hist );
364 twoF_pdf_hist = NULL;
369 while ( MC_iter-- ) {
383 REAL8 mism_rho2 = 0.0;
386 REAL8 prob_mismatch = 1.0;
388 REAL8 mismatch_from_pdf = 0.0;
389 REAL8 twoF_from_pdf = 0.0;
394 twoF = gsl_ran_flat( rng, 0, twoFs );
397 if ( calc_ABC_coeffs )
401 detector_states, noise_weights,
402 &A_coeff, &B_coeff, &C_coeff );
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 );
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 );
419 mismatch = gsl_ran_flat( rng, 0.0, max_mismatch );
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 );
427 prob_mismatch = prob_bin;
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;
440 if ( gsl_isnan( J ) || gsl_isnan( dJ ) ) {
447 REAL8 mism_pdf_cumul_prob_bin = gsl_ran_flat( rng, 0.0, 1.0 );
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;
458 mism_pdf_cumul_prob_bin -= prob_bin * ( right_bin - left_bin );
462 mism_rho2 = ( 1.0 - mismatch_from_pdf ) *
rho2;
468 if ( twoF_from_pdf < twoFs ) {
476 const size_t bin = twoF_from_pdf / twoF_pdf_hist_binw;
479 if ( !twoF_pdf_hist || bin >= twoF_pdf_hist->size )
480 if ( NULL == ( twoF_pdf_hist =
resize_histogram( twoF_pdf_hist, bin + 1 ) ) ) {
486 gsl_vector_int_set( twoF_pdf_hist, bin,
487 gsl_vector_int_get( twoF_pdf_hist, bin ) + 1 );
494 if ( gsl_isnan( J ) || gsl_isnan( dJ ) ) {
496 LogPrintf(
LOG_DEBUG,
"Reducing h0 to %0.4e and starting again because J=%0.4e, dJ=%0.4e\n",
h0, J, dJ );
498 MC_trials = MC_trials_init;
504 REAL8 MC_norm = twoFs / MC_trials;
506 MC_norm *= max_mismatch;
514 twoF_pdf_FDR = ( 1.0 * twoF_pdf_FD ) / MC_trials;
517 dh0 = ( FDR - J ) / dJ;
520 dh0 = GSL_SIGN( dh0 ) * GSL_MIN( fabs( dh0 ), fabs(
h0 * h0_brake ) );
523 fprintf(
fp,
" h0=%0.5e FDR_MC_int=%0.4f FDR_2F_dist=%0.4f\n",
h0, J, twoF_pdf_FDR );
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 );
528 if ( H0_ERROR <= max_rel_err || MC_trials >= MC_trials_reset ) {
538 MC_trials = (
INT8 )ceil( MC_trials * MC_trials_incr );
550 if ( MC_trials >= MC_trials_reset ) {
554 }
while ( MC_trials >= MC_trials_reset );
562 if ( ( fpH = fopen( twoF_pdf_hist_file,
"wb" ) ) == NULL ) {
563 XLALPrintError(
"Couldn't open histogram file '%s'\n", twoF_pdf_hist_file );
569 fprintf( fpH,
"%%%% %s\n%%%% %s\n",
"$Id$", cmdline );
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 ) );
587 gsl_matrix_free( mism_hist );
592 if ( twoF_pdf_hist ) {
593 gsl_vector_int_free( twoF_pdf_hist );
620 sky.longitude = gsl_ran_flat( rng,
alpha,
alpha + alpha_band );
621 sky.latitude = gsl_ran_flat( rng,
delta,
delta + delta_band );
638 return ( alpha_band != 0.0 || delta_band != 0.0 );
647 const REAL8 z = sqrt(
x * lambda );
650 gsl_error_handler_t *h = NULL;
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 );
658 gsl_set_error_handler( h );
661 return 0.5 * exp( -0.5 * (
x + lambda ) ) * sqrt(
x / lambda ) * I1.val;
670 const REAL8 z = sqrt(
x * lambda );
674 gsl_error_handler_t *h = NULL;
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 );
683 gsl_set_error_handler( h );
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 );
695 const REAL8 a = sqrt( lambda / 4 );
700 for (
i = 0;
i < 4; ++
i ) {
701 x += pow( gsl_ran_gaussian( rng, 1.0 ) +
a, 2.0 );
711 gsl_vector_int *new_hist = gsl_vector_int_alloc(
size );
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 );
int main(int argc, char *argv[])
REAL8 ran_ncx2_4(const gsl_rng *, REAL8)
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)
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)
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
static double double delta
const FstatOptionalArgs FstatOptionalArgsDefaults
Global initializer for setting FstatOptionalArgs to default values.
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.
MultiAMCoeffs * XLALComputeMultiAMCoeffs(const MultiDetectorStateSeries *multiDetStates, const MultiNoiseWeights *multiWeights, SkyPosition skypos)
Multi-IFO version of XLALComputeAMCoeffs().
#define XLAL_INIT_DECL(var,...)
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
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.
MultiNoiseWeights * XLALComputeMultiNoiseWeights(const MultiPSDVector *rngmed, UINT4 blocksRngMed, UINT4 excludePercentile)
Computes weight factors arising from MultiSFTs with different noise floors.
void XLALDestroyMultiNoiseWeights(MultiNoiseWeights *weights)
Destroy a MultiNoiseWeights object.
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
MultiSFTVector * XLALLoadMultiSFTs(const SFTCatalog *inputCatalog, REAL8 fMin, REAL8 fMax)
Function to load a catalog of SFTs from possibly different detectors.
void XLALDestroyMultiSFTVector(MultiSFTVector *multvect)
Destroy a multi SFT-vector.
SFTCatalog * XLALSFTdataFind(const CHAR *file_pattern, const SFTConstraints *constraints)
Find the list of SFTs matching the file_pattern and satisfying the given constraints,...
COORDINATESYSTEM_EQUATORIAL
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_MAIN(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK_NULL(assertion,...)
REAL8 Sinv_Tsft
normalization-factor (using single-sided PSD!)
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.
Multi-IFO container for antenna-pattern coefficients and atenna-pattern matrix .
AntennaPatternMatrix Mmunu
antenna-pattern matrix
Multi-IFO time-series of DetectorStates.
One noise-weight (number) per SFT (therefore indexed over IFOs and SFTs.
A collection of PSD vectors – one for each IFO in a multi-IFO search.
A collection of SFT vectors – one for each IFO in a multi-IFO search.
SFTVector ** data
sftvector for each ifo
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
SFTDescriptor * data
array of data-entries describing matched SFTs
UINT4 length
number of SFTs in catalog
'Constraints' for SFT-matching: which detector, within which time-stretch and which timestamps exactl...
SFTtype header
SFT-header info.