Commit 3e427a4e authored by Reinhard Prix's avatar Reinhard Prix
Browse files

Merge branch 'cfs_on_the_fly_make_noise2' into 'master'

Added support for on-the-fly noise generation in CFSv2

See merge request !178
parents bf8961b9 c26a76da
Pipeline #14539 passed with stages
in 22 minutes and 18 seconds
......@@ -182,6 +182,8 @@ typedef struct {
REAL8 dFreq;
CHAR transientOutputTimeUnit; /**< output format for transient times t0,tau: 'd' for days (deprecated) or 's' for seconds */
PulsarParamsVector *injectionSources; /**< Source parameters to inject: comma-separated list of file-patterns and/or direct config-strings ('{...}') */
MultiLIGOTimeGPSVector *multiTimestamps; /**< a vector of timestamps (only set if provided from dedicated time stamp files) */
MultiLALDetector multiIFO; /**< detectors to generate data for (if provided by user and not via noise files) */
} ConfigVariables;
......@@ -314,6 +316,11 @@ typedef struct {
BOOLEAN resampFFTPowerOf2; //!< in Resamp: enforce FFT length to be a power of two (by rounding up)
LALStringVector *injectionSources; /**< Source parameters to inject: comma-separated list of file-patterns and/or direct config-strings ('{...}') */
LALStringVector *injectSqrtSX; /**< Add Gaussian noise: list of respective detectors' noise-floors sqrt{Sn}" */
LALStringVector *IFOs; /**< list of detector-names "H1,H2,L1,.." or single detector */
LALStringVector *timestampsFiles; /**< Names of numDet timestamps files */
REAL8 Tsft; /**< length of one SFT in seconds, used in combination with timestamps files (otherwise taken from SFT files) */
INT4 randSeed; /**< allow user to specify random-number seed for reproducible noise-realizations */
// ----- deprecated and obsolete variables, kept around for backwards-compatibility -----
LIGOTimeGPS internalRefTime; /**< [DEPRECATED] internal reference time. Has no effect, XLALComputeFstat() now always uses midtime anyway ... */
......@@ -1010,6 +1017,11 @@ initUserVars ( UserInput_t *uvar )
uvar->transient_useFReg = 0;
uvar->resampFFTPowerOf2 = TRUE;
uvar->injectionSources = NULL;
uvar->injectSqrtSX = NULL;
uvar->IFOs = NULL;
uvar->timestampsFiles = NULL;
uvar->Tsft=1800.0;
/* ---------- register all user-variables ---------- */
XLALRegisterUvarMember( Alpha, RAJ, 'a', OPTIONAL, "Sky: equatorial J2000 right ascension (in radians or hours:minutes:seconds)");
......@@ -1053,7 +1065,7 @@ initUserVars ( UserInput_t *uvar )
XLALRegisterUvarMember( dorbitArgp, REAL8, 0, DEVELOPER, "Binary Orbit: Spacing in Orbital argument of periapse in radians");
XLALRegisterUvarMember( dorbitEcc, REAL8, 0, DEVELOPER, "Binary Orbit: Spacing in Orbital eccentricity");
XLALRegisterUvarMember(DataFiles, STRING, 'D', REQUIRED, "File-pattern specifying (also multi-IFO) input SFT-files");
XLALRegisterUvarMember(DataFiles, STRING, 'D', OPTIONAL, "File-pattern specifying (also multi-IFO) input SFT-files");
XLALRegisterUvarMember(IFO, STRING, 'I', OPTIONAL, "Detector: 'G1', 'L1', 'H1', 'H2' ...(useful for single-IFO v1-SFTs only!)");
XLALRegisterUvarMember( assumeSqrtSX, STRINGVector, 0, OPTIONAL, "Don't estimate noise-floors but assume (stationary) per-IFO sqrt{SX} (if single value: use for all IFOs)");
......@@ -1139,6 +1151,17 @@ initUserVars ( UserInput_t *uvar )
/* inject signals into the data being analyzed */
XLALRegisterUvarMember(injectionSources, STRINGVector, 0, DEVELOPER, "CSV list of files containing signal parameters for injection [see mfdv5]");
XLALRegisterUvarMember(injectSqrtSX, STRINGVector, 0, DEVELOPER, "Generate Gaussian Noise SFTs on-the-fly: CSV list of detectors' noise-floors sqrt{Sn}");
XLALRegisterUvarMember(IFOs, STRINGVector, 0, DEVELOPER, "CSV list of detectors, eg. \"H1,H2,L1,G1, ...\", when no SFT files are specified");
XLALRegisterUvarMember(timestampsFiles, STRINGVector, 0, DEVELOPER,
"Files containing timestamps for the generated SFTs when using "UVAR_STR( injectSqrtSX ) ". "
"Arguments correspond to the detectors given by " UVAR_STR( IFOs )
"; for example, if " UVAR_STR( IFOs ) " is set to 'H1,L1', then an argument of "
"'t1.txt,t2.txt' to this option will read H1 timestamps from the file 't1.txt', and L1 timestamps from the file 't2.txt'. "
"The timebase of the generated SFTs is specified by " UVAR_STR( Tsft ) ". "
);
XLALRegisterUvarMember(Tsft, REAL8, 0, DEVELOPER, "Generate SFTs with this timebase (in seconds) instead of loading from files. Requires --injectSqrtSX, --IFOs, --timestampsFiles");
XLALRegisterUvarMember(randSeed, INT4, 0, DEVELOPER, "Specify random-number seed for reproducible noise (0 means use /dev/urandom for seeding).");
// ---------- deprecated but still-supported or tolerated options ----------
XLALRegisterUvarMember ( RA, STRING, 0, DEPRECATED, "Use --Alpha instead" );
......@@ -1173,19 +1196,54 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar )
/* ----- set computational parameters for F-statistic from User-input ----- */
cfg->useResamp = ( uvar->FstatMethod >= FMETHOD_RESAMP_GENERIC ); // use resampling;
/* use IFO-contraint if one given by the user */
/* if IFO string vector was passed by user, parse it for later use */
if ( uvar->IFOs != NULL ) {
XLAL_CHECK ( XLALParseMultiLALDetector ( &(cfg->multiIFO), uvar->IFOs ) == XLAL_SUCCESS, XLAL_EFUNC );
}
/* read timestamps if requested by the user */
if (uvar->timestampsFiles != NULL) {
XLAL_CHECK ( (cfg->multiTimestamps = XLALReadMultiTimestampsFiles ( uvar->timestampsFiles )) != NULL, XLAL_EFUNC );
XLAL_CHECK ( (cfg->multiTimestamps->length > 0) && (cfg->multiTimestamps->data != NULL), XLAL_EINVAL, "Got empty timestamps-list from XLALReadMultiTimestampsFiles()\n" );
for ( UINT4 X=0; X < cfg->multiTimestamps->length; X ++ ) {
cfg->multiTimestamps->data[X]->deltaT = uvar->Tsft; // Tsft information not given by timestamps-file
}
}
/* use IFO-constraint if one given by the user */
if ( XLALUserVarWasSet ( &uvar->IFO ) ) {
XLAL_CHECK ( (constraints.detector = XLALGetChannelPrefix ( uvar->IFO )) != NULL, XLAL_EFUNC );
}
LIGOTimeGPS minStartTime = uvar->minStartTime;
LIGOTimeGPS maxStartTime = uvar->maxStartTime;
constraints.minStartTime = &minStartTime;
constraints.maxStartTime = &maxStartTime;
/* get full SFT-catalog of all matching (multi-IFO) SFTs */
LogPrintf (LOG_NORMAL, "Finding all SFTs to load ... ");
XLAL_CHECK ( (catalog = XLALSFTdataFind ( uvar->DataFiles, &constraints )) != NULL, XLAL_EFUNC );
LogPrintfVerbatim (LOG_NORMAL, "done. (found %d SFTs)\n", catalog->length);
/* DataFiles optional because of injectSqrtSX option, don't try to load if no files given **/
if( uvar->DataFiles != NULL ) {
LogPrintf (LOG_NORMAL, "Finding all SFTs to load ... ");
XLAL_CHECK ( (catalog = XLALSFTdataFind ( uvar->DataFiles, &constraints )) != NULL, XLAL_EFUNC );
LogPrintfVerbatim (LOG_NORMAL, "done. (found %d SFTs)\n", catalog->length);
} else {
/* Build a fake catalog with timestamps and IFOs given on the commandline instead of noise data files */
/* the data missing in the locators then signal to the Fstat code that fake noise needs to be generated, */
/* the noise level is passed from the sqrtSX option */
XLAL_CHECK ( (catalog = XLALMultiAddToFakeSFTCatalog ( NULL, uvar->IFOs,cfg-> multiTimestamps)) != NULL, XLAL_EFUNC );
}
if ( constraints.detector ) {
XLALFree ( constraints.detector );
......@@ -1195,6 +1253,7 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar )
/* deduce start- and end-time of the observation spanned by the data */
UINT4 numSFTfiles = catalog->length;
cfg->Tsft = 1.0 / catalog->data[0].header.deltaF;
cfg->startTime = catalog->data[0].header.epoch;
endTime = catalog->data[numSFTfiles - 1].header.epoch;
......@@ -1415,17 +1474,26 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar )
}
MultiNoiseFloor *injectSqrtSX = NULL;
MultiNoiseFloor s_injectSqrtSX;
if(uvar->injectSqrtSX != NULL) {
XLAL_CHECK( XLALParseMultiNoiseFloor( &s_injectSqrtSX, uvar->injectSqrtSX, XLALCountIFOsInCatalog(catalog) ) == XLAL_SUCCESS, XLAL_EFUNC );
injectSqrtSX = &s_injectSqrtSX;
}
FstatOptionalArgs optionalArgs = FstatOptionalArgsDefaults;
optionalArgs.Dterms = uvar->Dterms;
optionalArgs.SSBprec = uvar->SSBprecision;
optionalArgs.runningMedianWindow = uvar->RngMedWindow;
optionalArgs.injectSources = cfg->injectionSources;
optionalArgs.injectSqrtSX = injectSqrtSX;
optionalArgs.randSeed=uvar->randSeed;
optionalArgs.assumeSqrtSX = assumeSqrtSX;
optionalArgs.FstatMethod = uvar->FstatMethod;
optionalArgs.resampFFTPowerOf2 = uvar->resampFFTPowerOf2;
optionalArgs.collectTiming = XLALUserVarWasSet ( &uvar->outputFstatTiming );
XLAL_CHECK ( (cfg->Fstat_in = XLALCreateFstatInput( catalog, fCoverMin, fCoverMax, cfg->dFreq, cfg->ephemeris, &optionalArgs )) != NULL, XLAL_EFUNC );
XLALDestroySFTCatalog(catalog);
......@@ -1744,6 +1812,10 @@ Freemem( ConfigVariables *cfg )
XLALDestroyPulsarParamsVector ( cfg->injectionSources );
}
if( cfg->multiTimestamps ) {
XLALDestroyMultiTimestamps( cfg->multiTimestamps );
}
return;
} /* Freemem() */
......@@ -1935,6 +2007,22 @@ checkUserInputConsistency ( const UserInput_t *uvar )
/* check SignalOnly and assumeSqrtSX */
XLAL_CHECK ( !uvar->SignalOnly || (uvar->assumeSqrtSX == NULL), XLAL_EINVAL, "Cannot pass --SignalOnly AND --assumeSqrtSX at the same time!\n");
XLAL_CHECK ( (uvar->DataFiles ==NULL) ^ (uvar->injectSqrtSX == NULL), XLAL_EINVAL, "Must pass exactly one out of --DataFiles or --injectSqrtSX \n");
if(uvar->DataFiles !=NULL) {
XLAL_CHECK ( !XLALUserVarWasSet(&uvar->Tsft) , XLAL_EINVAL, UVAR_STR(Tsft) " can only be used for data generation with " UVAR_STR(injectSqrtSX)", not when loading existing " UVAR_STR(DataFiles) "\n");
XLAL_CHECK ( uvar->IFOs == NULL , XLAL_EINVAL, UVAR_STR(IFOs) " can only be used for data generation with " UVAR_STR(injectSqrtSX) ", not when loading existing " UVAR_STR(DataFiles) "\n");
XLAL_CHECK ( uvar->timestampsFiles == NULL , XLAL_EINVAL,UVAR_STR(timestampsFiles) " can only be used for data generation with " UVAR_STR(injectSqrtSX) ", not when loading existing " UVAR_STR(DataFiles) "\n");
}
if(uvar->injectSqrtSX !=NULL) {
XLAL_CHECK ( uvar->timestampsFiles != NULL && uvar->IFOs != NULL , XLAL_EINVAL,"--injectSqrtSX requires --IFOs, --timestampsFiles \n");
}
return XLAL_SUCCESS;
} /* checkUserInputConsistency() */
......
%% LAL: 6.18.0.1 (UNCLEAN a9590f59eb2e6a58cc9087a496eda4cc79cfead0)
%% LALSimulation: 1.7.3.1 (UNCLEAN a9590f59eb2e6a58cc9087a496eda4cc79cfead0)
%% LALPulsar: 1.16.0.1 (UNCLEAN a9590f59eb2e6a58cc9087a496eda4cc79cfead0)
%% LALApps: 6.21.0.1 (UNCLEAN a9590f59eb2e6a58cc9087a496eda4cc79cfead0)
%% cmdline: lalapps_dumpSFT --SFTfiles="H-80_H1_1800SFT_mfdv5-711595934-144000.sft" --timestampsOnly=TRUE
%% Timestamps: seconds nano-seconds
711595934 000000000
711597734 000000000
711599534 000000000
711601334 000000000
711603134 000000000
711604934 000000000
711606734 000000000
711608534 000000000
711610334 000000000
711612134 000000000
711613934 000000000
711615734 000000000
711617534 000000000
711619334 000000000
711621134 000000000
711622934 000000000
711624734 000000000
711626534 000000000
711628334 000000000
711630134 000000000
711631934 000000000
711633734 000000000
711635534 000000000
711637334 000000000
711639134 000000000
711640934 000000000
711642734 000000000
711644534 000000000
711646334 000000000
711648134 000000000
711649934 000000000
711651734 000000000
711653534 000000000
711655334 000000000
711657134 000000000
711658934 000000000
711660734 000000000
711662534 000000000
711664334 000000000
711666134 000000000
711667934 000000000
711669734 000000000
711671534 000000000
711673334 000000000
711675134 000000000
711676934 000000000
711678734 000000000
711680534 000000000
711682334 000000000
711684134 000000000
711685934 000000000
711687734 000000000
711689534 000000000
711691334 000000000
711693134 000000000
711694934 000000000
711696734 000000000
711698534 000000000
711700334 000000000
711702134 000000000
711703934 000000000
711705734 000000000
711707534 000000000
711709334 000000000
711711134 000000000
711712934 000000000
711714734 000000000
711716534 000000000
711718334 000000000
711720134 000000000
711721934 000000000
711723734 000000000
711725534 000000000
711727334 000000000
711729134 000000000
711730934 000000000
711732734 000000000
711734534 000000000
711736334 000000000
711738134 000000000
%% LAL: 6.18.0.1 (UNCLEAN a9590f59eb2e6a58cc9087a496eda4cc79cfead0)
%% LALSimulation: 1.7.3.1 (UNCLEAN a9590f59eb2e6a58cc9087a496eda4cc79cfead0)
%% LALPulsar: 1.16.0.1 (UNCLEAN a9590f59eb2e6a58cc9087a496eda4cc79cfead0)
%% LALApps: 6.21.0.1 (UNCLEAN a9590f59eb2e6a58cc9087a496eda4cc79cfead0)
%% cmdline: lalapps_dumpSFT --SFTfiles="L-80_L1_1800SFT_mfdv5-711595934-144000.sft" --timestampsOnly=TRUE
%% Timestamps: seconds nano-seconds
711595934 000000000
711597734 000000000
711599534 000000000
711601334 000000000
711603134 000000000
711604934 000000000
711606734 000000000
711608534 000000000
711610334 000000000
711612134 000000000
711613934 000000000
711615734 000000000
711617534 000000000
711619334 000000000
711621134 000000000
711622934 000000000
711624734 000000000
711626534 000000000
711628334 000000000
711630134 000000000
711631934 000000000
711633734 000000000
711635534 000000000
711637334 000000000
711639134 000000000
711640934 000000000
711642734 000000000
711644534 000000000
711646334 000000000
711648134 000000000
711649934 000000000
711651734 000000000
711653534 000000000
711655334 000000000
711657134 000000000
711658934 000000000
711660734 000000000
711662534 000000000
711664334 000000000
711666134 000000000
711667934 000000000
711669734 000000000
711671534 000000000
711673334 000000000
711675134 000000000
711676934 000000000
711678734 000000000
711680534 000000000
711682334 000000000
711684134 000000000
711685934 000000000
711687734 000000000
711689534 000000000
711691334 000000000
711693134 000000000
711694934 000000000
711696734 000000000
711698534 000000000
711700334 000000000
711702134 000000000
711703934 000000000
711705734 000000000
711707534 000000000
711709334 000000000
711711134 000000000
711712934 000000000
711714734 000000000
711716534 000000000
711718334 000000000
711720134 000000000
711721934 000000000
711723734 000000000
711725534 000000000
711727334 000000000
711729134 000000000
711730934 000000000
711732734 000000000
711734534 000000000
711736334 000000000
711738134 000000000
......@@ -78,7 +78,9 @@ EXTRA_DIST += \
testCFSv2_grid3.dat.ref.gz \
testCFSv2_grid6.dat.ref.gz \
testCFSv2_grid8.dat.ref.gz \
testCFSv2_grid9.dat.ref.gz
testCFSv2_grid9.dat.ref.gz \
H1_test_timestamps.dat \
L1_test_timestamps.dat
TESTS = testPredictFstat.sh \
testSynthLV.sh
......
......@@ -93,9 +93,15 @@ loudest_Demod=${testDir}/Loudest_Demod.dat
outfile_Resamp=${testDir}/Fstat_Resamp.dat
loudest_Resamp=${testDir}/Loudest_Resamp.dat
outfile_Resamp_otfn=${testDir}/Fstat_Resamp_otfn.dat
loudest_Resamp_otfn=${testDir}/Loudest_Resamp_otfn.dat
timefile_Comp=${testDir}/timing_Comp.dat
timefile_Demod=${testDir}/timing_Demod.dat
timefile_Resamp=${testDir}/timing_Resamp.dat
timefile_Resamp_otfn=${testDir}/timing_Resamp_otfn.dat
##--------------------------------------------------
## test starts here
......@@ -165,13 +171,31 @@ if ! eval "$cmdline 2> /dev/null"; then
exit 1;
fi
echo
echo "----------------------------------------------------------------------"
echo " STEP 5: run directed CFS_v2 with resampling and on-the-fly noise"
echo "----------------------------------------------------------------------"
echo
cmdline="$cfs_code $cfs_CL --injectionSources='${injectionSources}' --injectSqrtSX=${sqrtSX} --IFOs='H1,L1' --Tsft=${Tsft} --randSeed=1 --timestampsFiles='H1_test_timestamps.dat,L1_test_timestamps.dat' --FstatMethod=ResampBest --outputFstat=$outfile_Resamp_otfn --outputTiming=$timefile_Resamp_otfn --outputLoudest=${loudest_Resamp_otfn}"
echo $cmdline;
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$cfs_code' ..."
exit 1;
fi
echo "----------------------------------------"
echo " STEP 5: Comparing results: "
echo " STEP 6: Comparing results: "
echo "----------------------------------------"
sort -o $outfile_Comp $outfile_Comp
sort -o $outfile_Demod $outfile_Demod
sort -o $outfile_Resamp $outfile_Resamp
sort -o $outfile_Resamp_otfn $outfile_Resamp_otfn
## compare absolute differences instead of relative, allow deviations of up to sigma=sqrt(8)~2.8
echo
......@@ -194,9 +218,23 @@ else
echo " ==> OK."
fi
#laxer tolerances as we have different noiuse realizations in this testcase. Incuded just for smoke test
echo
cmdline="$cmp_code -1 ./${outfile_Resamp} -2 ./${outfile_Resamp_otfn} --tol-L1=7e-1 --tol-L2=7e-1 --tol-angle=0.2 --tol-atMax=2e-1"
echo -n $cmdline
if ! eval $cmdline; then
echo "==> OUCH... files differ. Something might be wrong with --injectSqrtSX et al."
exit 2
else
echo " ==> OK."
fi
echo
echo "----------------------------------------------------------------------"
echo " STEP 6: Sanity-check Resampling parameter estimation: "
echo " STEP 7: Sanity-check Resampling parameter estimation: "
echo "----------------------------------------------------------------------"
echo
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment