Commit b2facb39 authored by Vivien Raymond's avatar Vivien Raymond
Browse files

Merge branch 'master' into cbc_bayesian_devel

Original: 1f6b21e7a92fa415333258fb89a2f40e63bcf259
parents ac91d690 bca74b0f
......@@ -18,6 +18,7 @@ AC_SUBST([abs_lalsuite_srcdir],[`cd ${srcdir} && pwd`])
AC_PROG_LN_S
AC_CHECK_PROGS([GIT], [git], [])
LALSUITE_ENABLE_ALL_LAL
LALSUITE_ENABLE_LALFRAME
LALSUITE_ENABLE_LALMETAIO
LALSUITE_ENABLE_LALXML
......
# lalsuite_build.m4 - top level build macros
#
# serial 15
# serial 16
AC_DEFUN([LALSUITE_USE_LIBTOOL],
[## $0: Generate a libtool script for use in configure tests
......@@ -129,6 +129,18 @@ AC_DEFUN([LALSUITE_ENABLE_NIGHTLY],
AC_SUBST(NIGHTLY_VERSION)
])
AC_DEFUN([LALSUITE_ENABLE_ALL_LAL],
[AC_ARG_ENABLE(
[all_lal],
AC_HELP_STRING([--enable-all-lal],[enable/disable compilation of all LAL libraries]),
[ case "${enableval}" in
yes) all_lal=true;;
no) all_lal=false;;
*) AC_MSG_ERROR(bad value ${enableval} for --enable-all-lal) ;;
esac
], [ all_lal= ] )
])
AC_DEFUN([LALSUITE_ENABLE_LALFRAME],
[AC_ARG_ENABLE(
[lalframe],
......@@ -138,7 +150,7 @@ AC_DEFUN([LALSUITE_ENABLE_LALFRAME],
no) lalframe=false;;
*) AC_MSG_ERROR(bad value ${enableval} for --enable-lalframe) ;;
esac
], [ lalframe=true ] )
], [ lalframe=${all_lal:-true} ] )
if test "$frame" = "false"; then
lalframe=false
fi
......@@ -153,7 +165,7 @@ AC_DEFUN([LALSUITE_ENABLE_LALMETAIO],
no) lalmetaio=false;;
*) AC_MSG_ERROR(bad value ${enableval} for --enable-lalmetaio) ;;
esac
], [ lalmetaio=true ] )
], [ lalmetaio=${all_lal:-true} ] )
if test "$metaio" = "false"; then
lalmetaio=false
fi
......@@ -168,7 +180,7 @@ AC_DEFUN([LALSUITE_ENABLE_LALXML],
no) lalxml=false;;
*) AC_MSG_ERROR(bad value ${enableval} for --enable-lalxml) ;;
esac
], [ lalxml=false ] )
], [ lalxml=${all_lal:-false} ] )
])
AC_DEFUN([LALSUITE_ENABLE_LALBURST],
......@@ -180,7 +192,7 @@ AC_DEFUN([LALSUITE_ENABLE_LALBURST],
no) lalburst=false;;
*) AC_MSG_ERROR(bad value ${enableval} for --enable-lalburst) ;;
esac
], [ lalburst=true ] )
], [ lalburst=${all_lal:-true} ] )
if test "$lalmetaio" = "false"; then
lalburst=false
fi])
......@@ -194,7 +206,7 @@ AC_DEFUN([LALSUITE_ENABLE_LALINSPIRAL],
no) lalinspiral=false;;
*) AC_MSG_ERROR(bad value ${enableval} for --enable-lalinspiral) ;;
esac
], [ lalinspiral=true ] )
], [ lalinspiral=${all_lal:-true} ] )
if test "$lalmetaio" = "false"; then
lalinspiral=false
fi
......@@ -209,7 +221,7 @@ AC_DEFUN([LALSUITE_ENABLE_LALPULSAR],
no) lalpulsar=false;;
*) AC_MSG_ERROR(bad value ${enableval} for --enable-lalpulsar) ;;
esac
], [ lalpulsar=true ] )
], [ lalpulsar=${all_lal:-true} ] )
])
AC_DEFUN([LALSUITE_ENABLE_LALSTOCHASTIC],
......@@ -221,7 +233,7 @@ AC_DEFUN([LALSUITE_ENABLE_LALSTOCHASTIC],
no) lalstochastic=false;;
*) AC_MSG_ERROR(bad value ${enableval} for --enable-lalstochastic) ;;
esac
], [ lalstochastic=true ] )
], [ lalstochastic=${all_lal:-true} ] )
if test "$lalmetaio" = "false"; then
lalstochastic=false
fi
......
......@@ -123,7 +123,7 @@ int XLALComputeDQ(REAL4* sv_data, int r_sv,
badgamma = 0;
for (j = 0; j < r_gamma; j++) {
REAL8 re = gamma_data[j].re;
REAL8 re = gamma_data[i*r_gamma + j].re;
if (re < 0.8 || re > 1.2) /* || isnan(re) || isinf(re) not C89 */
badgamma = 1;
}
......
......@@ -44,6 +44,9 @@
#include <lal/NormalizeSFTRngMed.h>
#include <lal/LALInitBarycenter.h>
#include <lal/SFTClean.h>
#include <lal/Segments.h>
#include <lal/FrequencySeries.h>
#include <lal/Units.h>
#include <lal/LogPrintf.h>
......@@ -103,8 +106,8 @@ typedef struct
REAL8 Freq; /**< *physical* start frequency to compute PSD for (excluding rngmed wings) */
REAL8 FreqBand; /**< *physical* frequency band to compute PSD for (excluding rngmed wings) */
REAL8 startTime;
REAL8 endTime;
REAL8 startTime; /**< earliest SFT-timestamp to include */
REAL8 endTime; /**< last SFT-timestamp to include */
CHAR *IFO;
CHAR *timeStampsFile;
LALStringVector *linefiles;
......@@ -126,15 +129,26 @@ typedef struct
BOOLEAN outFreqBinEnd; /**< output the end frequency of each bin? */
BOOLEAN dumpMultiPSDVector; /**< output multi-PSD vector over IFOs, timestamps, and frequencies into file(s) */
CHAR *outputQ; /**< output the 'data-quality factor' Q(f) into this file */
REAL8 fStart; /**< Start Frequency to load from SFT and compute PSD, including wings (it is RECOMMENDED to use --Freq instead) */
REAL8 fBand; /**< Frequency Band to load from SFT and compute PSD, including wings (it is RECOMMENDED to use --FreqBand instead) */
} UserVariables_t;
/** Config variables 'derived' from user-input
*/
typedef struct
{
UINT4 firstBin; /**< first PSD bin for output */
UINT4 lastBin; /**< last PSD bin for output */
LALSeg dataSegment; /**< the data-segment for which PSD was computed */
} ConfigVariables_t;
/*---------- empty structs for initializations ----------*/
UserVariables_t empty_UserVariables;
ConfigVariables_t empty_ConfigVariables;
/* ---------- global variables ----------*/
extern int vrbflg;
......@@ -146,9 +160,11 @@ int initUserVars (int argc, char *argv[], UserVariables_t *uvar);
void LALfwriteSpectrograms ( LALStatus *status, const CHAR *bname, const MultiPSDVector *multiPSD );
static REAL8 math_op(REAL8*, size_t, INT4);
int XLALDumpMultiPSDVector ( const CHAR *outbname, const MultiPSDVector *multiPSDVect );
MultiSFTVector *XLALReadSFTs ( UINT4 *firstBin, UINT4 *lastBin, const UserVariables_t *uvar );
MultiSFTVector *XLALReadSFTs ( ConfigVariables_t *cfg, const UserVariables_t *uvar );
int XLALCropMultiPSDVector ( MultiPSDVector *multiPSDVect, UINT4 firstBin, UINT4 lastBin );
REAL8FrequencySeries *XLALComputeSegmentDataQ ( const MultiPSDVector *multiPSDVect, LALSeg segment );
int XLALWriteREAL8FrequencySeries_to_file ( const REAL8FrequencySeries *series, const char *fname );
/*============================================================
* FUNCTION definitions
......@@ -159,6 +175,7 @@ main(int argc, char *argv[])
const char *fn = __func__;
static LALStatus status; /* LALStatus pointer */
UserVariables_t uvar = empty_UserVariables;
ConfigVariables_t cfg = empty_ConfigVariables;
UINT4 k, numBins, numIFOs, maxNumSFTs, X, alpha;
REAL8 Freq0, dFreq, normPSD;
......@@ -189,9 +206,8 @@ main(int argc, char *argv[])
if (uvar.help)
return EXIT_SUCCESS;
UINT4 firstBin, lastBin;
MultiSFTVector *inputSFTs = NULL;
if ( ( inputSFTs = XLALReadSFTs ( &firstBin, &lastBin, &uvar ) ) == NULL )
if ( ( inputSFTs = XLALReadSFTs ( &cfg, &uvar ) ) == NULL )
{
XLALPrintError ("Call to XLALReadSFTs() failed with xlalErrno = %d\n", xlalErrno );
return EXIT_FAILURE;
......@@ -227,8 +243,8 @@ main(int argc, char *argv[])
MultiPSDVector *multiPSD = NULL;
LAL_CALL( LALNormalizeMultiSFTVect (&status, &multiPSD, inputSFTs, uvar.blocksRngMed), &status);
/* restrict this PSD to just the "physical" band if requested using {--Freq, --FreqBand} */
if ( ( XLALCropMultiPSDVector ( multiPSD, firstBin, lastBin )) != XLAL_SUCCESS ) {
XLALPrintError ("%s: XLALCropMultiPSDVector (inputPSD, %d, %d) failed with xlalErrno = %d\n", fn, firstBin, lastBin, xlalErrno );
if ( ( XLALCropMultiPSDVector ( multiPSD, cfg.firstBin, cfg.lastBin )) != XLAL_SUCCESS ) {
XLALPrintError ("%s: XLALCropMultiPSDVector (inputPSD, %d, %d) failed with xlalErrno = %d\n", fn, cfg.firstBin, cfg.lastBin, xlalErrno );
return EXIT_FAILURE;
}
......@@ -346,6 +362,19 @@ main(int argc, char *argv[])
}
} /* if uvar.dumpMultiPSDVector */
/* ----- if requested, compute data-quality factor 'Q' -------------------- */
if ( uvar.outputQ )
{
REAL8FrequencySeries *Q;
if ( (Q = XLALComputeSegmentDataQ ( multiPSD, cfg.dataSegment )) == NULL ) {
XLALPrintError ("%s: XLALComputeSegmentDataQ() failed with xlalErrno = %d\n", fn, xlalErrno );
return EXIT_FAILURE;
}
if ( XLAL_SUCCESS != XLALWriteREAL8FrequencySeries_to_file ( Q, uvar.outputQ ) ) {
return EXIT_FAILURE;
}
XLALDestroyREAL8FrequencySeries ( Q );
} /* if outputQ */
/* ---------- BINNING if requested ---------- */
/* work out bin size */
......@@ -559,13 +588,14 @@ initUserVars (int argc, char *argv[], UserVariables_t *uvar)
XLALregBOOLUserStruct (help, 'h', UVAR_HELP, "Print this message" );
XLALregSTRINGUserStruct(inputData, 'i', UVAR_REQUIRED, "Input SFT pattern");
XLALregSTRINGUserStruct(outputPSD, 'o', UVAR_OPTIONAL, "Output PSD into this file");
XLALregSTRINGUserStruct(outputQ, 0, UVAR_OPTIONAL, "Output the 'data-quality factor' Q(f) into this file");
XLALregSTRINGUserStruct(outputSpectBname, 0 , UVAR_OPTIONAL, "Filename-base for (binary) spectrograms (one per IFO)");
XLALregREALUserStruct (Freq, 0, UVAR_OPTIONAL, "physical start frequency to compute PSD for (excluding rngmed wings)");
XLALregREALUserStruct (FreqBand, 0, UVAR_OPTIONAL, "physical frequency band to compute PSD for (excluding rngmed wings)");
XLALregREALUserStruct (startTime, 's', UVAR_OPTIONAL, "GPS start time");
XLALregREALUserStruct (endTime, 'e', UVAR_OPTIONAL, "GPS end time");
XLALregREALUserStruct (startTime, 's', UVAR_OPTIONAL, "GPS timestamp of earliest SFT to include");
XLALregREALUserStruct (endTime, 'e', UVAR_OPTIONAL, "GPS timestamp of last SFT to include (NOTE: this refers to the SFT start-time!)");
XLALregSTRINGUserStruct(timeStampsFile, 't', UVAR_OPTIONAL, "Time-stamps file");
XLALregSTRINGUserStruct(IFO, 0 , UVAR_OPTIONAL, "Detector filter");
......@@ -888,15 +918,16 @@ XLALDumpMultiPSDVector ( const CHAR *outbname, /**< output basename <outbname>
/** Load all SFTs according to user-input, returns multi-SFT vector.
*
* \return cfg:
* Returns 'effective' range of SFT-bins [firstBin, lastBin], which which the PSD will be estimated:
* - if the user input {fStart, fBand} then these are loaded from SFTs and directly translated into bins
* - if user input {Freq, FreqBand}, we load a wider frequency-band ADDING running-median/2 on either side
* from the SFTs, and firstBind, lastBin correspond to {Freq,FreqBand} (rounded to closest bins)
* Also returns the 'data-segment' for which SFTs were loaded
*
*/
MultiSFTVector *
XLALReadSFTs ( UINT4 *firstBin, /**< [out] first PSD bin for output */
UINT4 *lastBin, /**< [out] last PSD bin for output */
XLALReadSFTs ( ConfigVariables_t *cfg, /**< [out] return derived configuration info (firstBin, lastBin, segment) */
const UserVariables_t *uvar /**< [in] complete user-input */
)
{
......@@ -904,7 +935,7 @@ XLALReadSFTs ( UINT4 *firstBin, /**< [out] first PSD bin for output */
SFTCatalog *catalog = NULL;
SFTConstraints constraints = empty_SFTConstraints;
LIGOTimeGPS startTimeGPS, endTimeGPS;
LIGOTimeGPS startTimeGPS = {0,0}, endTimeGPS = {0,0};
LIGOTimeGPSVector *inputTimeStampsVector = NULL;
/* check input */
......@@ -912,6 +943,10 @@ XLALReadSFTs ( UINT4 *firstBin, /**< [out] first PSD bin for output */
XLALPrintError ("%s: invalid NULL input 'uvar' or 'uvar->inputData'\n", fn );
XLAL_ERROR_NULL ( fn, XLAL_EINVAL );
}
if ( !cfg ) {
XLALPrintError ("%s: invalid NULL input 'cfg'", fn );
XLAL_ERROR_NULL ( fn, XLAL_EINVAL );
}
/* set detector constraint */
if ( XLALUserVarWasSet ( &uvar->IFO ) )
......@@ -974,7 +1009,7 @@ XLALReadSFTs ( UINT4 *firstBin, /**< [out] first PSD bin for output */
/* ---------- figure out the right frequency-band to read from the SFTs, depending on user-input ----- */
REAL8 fMin, fMax;
UINT4 binsOffset; /* rngmed bin offset from start and end */
UINT4 binsBand; /* width of physical FreqBand in bins */
UINT4 binsBand=0; /* width of physical FreqBand in bins */
if ( have_Freq )
{
REAL8 dFreq = catalog->data[0].header.deltaF;
......@@ -993,6 +1028,16 @@ XLALReadSFTs ( UINT4 *firstBin, /**< [out] first PSD bin for output */
binsOffset = 0; /* no truncation of rngmed sidebands */
}
/* ----- figure out the data-segment span from the user-input and SFT-catalog ----- */
/* if used passed these, then 'startTimeGPS' and 'endTimeGPS' are already set */
if ( startTimeGPS.gpsSeconds == 0 )
startTimeGPS = catalog->data[0].header.epoch;
if ( endTimeGPS.gpsSeconds == 0 )
endTimeGPS = catalog->data[catalog->length-1].header.epoch;
/* SFT 'constraints' only refer to SFT *start-times*, for segment we need the end-time */
REAL8 Tsft = 1.0 / catalog->data[0].header.deltaF;
XLALGPSAdd ( &endTimeGPS, Tsft );
/* ---------- read the sfts ---------- */
LogPrintf (LOG_DEBUG, "Loading all SFTs ... ");
MultiSFTVector *multi_sfts;
......@@ -1019,8 +1064,10 @@ XLALReadSFTs ( UINT4 *firstBin, /**< [out] first PSD bin for output */
}
/* return results */
(*firstBin) = (UINT4) bin0;
(*lastBin ) = (UINT4) bin1;
cfg->firstBin = (UINT4) bin0;
cfg->lastBin = (UINT4) bin1;
cfg->dataSegment.start = startTimeGPS;
cfg->dataSegment.end = endTimeGPS;
XLALPrintInfo ("%s: loaded SFTs have %d bins, effective PSD output band is [%d, %d]\n", fn, numBins, bin0, bin1 );
......@@ -1107,3 +1154,175 @@ XLALCropMultiPSDVector ( MultiPSDVector *multiPSDVect,
return XLAL_SUCCESS;
} /* XLALCropMultiPSDVector() */
/** Compute the "data-quality factor" \f$\mathcal{Q}(f) = \sum_X \frac{\epsilon_X}{\mathcal{S}_X(f)}\f$ over the given SFTs.
* The input \a multiPSD is a pre-computed PSD map \f$\mathcal{S}_{X,i}(f)\f$, over IFOs \f$X\f$, SFTs \f$i\f$
* and frequency \f$f\f$.
*
* \return the output is a vector \f$\mathcal{Q}(f)\f$.
*
*/
REAL8FrequencySeries *
XLALComputeSegmentDataQ ( const MultiPSDVector *multiPSDVect, /**< input PSD map over IFOs, SFTs, and frequencies */
LALSeg segment /**< segment to compute Q for */
)
{
const char *fn = __func__;
/* check input consistency */
if ( multiPSDVect == NULL ) {
XLALPrintError ("%s: NULL input 'multiPSDVect'\n", fn );
XLAL_ERROR_NULL ( fn, XLAL_EINVAL );
}
if ( multiPSDVect->length == 0 || multiPSDVect->data==0 ) {
XLALPrintError ("%s: invalid multiPSDVect input (length=0 or data=NULL)\n", fn );
XLAL_ERROR_NULL ( fn, XLAL_EINVAL );
}
REAL8 Tseg = XLALGPSDiff ( &segment.end, &segment.start );
if ( Tseg <= 0 ) {
XLALPrintError ("%s: negative segment-duration '%g'\n", fn, Tseg );
XLAL_ERROR_NULL ( fn, XLAL_EINVAL );
}
REAL8 Tsft = 1.0 / multiPSDVect->data[0]->data[0].deltaF;
REAL8 f0 = multiPSDVect->data[0]->data[0].f0;
REAL8 dFreq = multiPSDVect->data[0]->data[0].deltaF;
UINT4 numFreqs = multiPSDVect->data[0]->data[0].data->length;
REAL8FrequencySeries *Q, *SXinv;
if ( (Q = XLALCreateREAL8FrequencySeries ( "Qfactor", &segment.start, f0, dFreq, &lalHertzUnit, numFreqs )) == NULL ) {
XLALPrintError ("%s: Q = XLALCreateREAL8FrequencySeries(numFreqs=%d) failed with xlalErrno = %d\n", fn, numFreqs, xlalErrno );
XLAL_ERROR_NULL ( fn, XLAL_EFUNC );
}
if ( (SXinv = XLALCreateREAL8FrequencySeries ( "SXinv", &segment.start, f0, dFreq, &lalHertzUnit, numFreqs )) == NULL ) {
XLALPrintError ("%s: SXinv = XLALCreateREAL8FrequencySeries(numFreqs=%d) failed with xlalErrno = %d\n", fn, numFreqs, xlalErrno );
XLAL_ERROR_NULL ( fn, XLAL_EFUNC );
}
/* initialize Q-array to zero, as we'll keep adding to it */
memset ( Q->data->data, 0, Q->data->length * sizeof(Q->data->data[0]) );
/* ----- loop over IFOs ----- */
UINT4 numIFOs = multiPSDVect->length;
UINT4 X;
for ( X = 0; X < numIFOs; X ++ )
{
PSDVector *thisPSDVect = multiPSDVect->data[X];
/* initialize SXinv-array to zero, as we'll keep adding to it */
memset ( SXinv->data->data, 0, SXinv->data->length * sizeof(SXinv->data->data[0]) );
UINT4 numSFTsInSeg = 0; /* reset counter of SFTs within this segment */
/* ----- loop over all timestamps ----- */
/* find SFTs inside segment, count them and combine their PSDs */
UINT4 numTS = thisPSDVect->length;
UINT4 iTS;
for ( iTS = 0; iTS < numTS; iTS++ )
{
REAL8FrequencySeries *thisPSD = &thisPSDVect->data[iTS];
/* some internal consistency/paranoia checks */
if ( ( f0 != thisPSD->f0) || ( dFreq != thisPSD->deltaF ) || (numFreqs != thisPSD->data->length ) ) {
XLALPrintError ("%s: %d-th timestamp %f: inconsistent PSDVector: f0 = %g : %g, dFreq = %g : %g, numFreqs = %d : %d \n",
fn, iTS, XLALGPSGetREAL8( &thisPSD->epoch ), f0, thisPSD->f0, dFreq, thisPSD->deltaF, numFreqs, thisPSD->data->length );
XLAL_ERROR_NULL ( fn, XLAL_EDOM );
}
LIGOTimeGPS gpsStart = thisPSD->epoch;
LIGOTimeGPS gpsEnd = gpsStart;
XLALGPSAdd( &gpsEnd, Tsft - 1e-3 ); /* subtract 1ms from end: segments are half-open intervals [t0, t1) */
int cmp1 = XLALGPSInSeg ( &gpsStart, &segment );
int cmp2 = XLALGPSInSeg ( &gpsEnd, &segment );
if ( cmp2 < 0 ) /* SFT-end before segment => advance to the next one */
continue;
if ( cmp1 > 0 ) /* SFT-start past end of segment: ==> terminate loop */
break;
if ( (cmp1 == 0) && (cmp2 == 0) ) /* this SFT is inside segment */
{
numSFTsInSeg ++;
/* add SXinv(f) += 1/SX_i(f) over all frequencies */
UINT4 iFreq;
for ( iFreq = 0; iFreq < numFreqs; iFreq++ )
SXinv->data->data[iFreq] += 1.0 / thisPSD->data->data[iFreq] ;
} /* if SFT inside segment */
} /* for iTS < numTS */
/* compute duty-cycle eps_X = nX * Tsft / Tseg for this IFO */
REAL8 duty_X = numSFTsInSeg * Tsft / Tseg;
/* sanity check: eps in [0, 1]*/
if ( (duty_X < 0) && (duty_X > 1 ) ) {
XLALPrintError ("%s: something is WRONG: duty-cyle = %g not within [0,1]!\n", fn, duty_X );
XLAL_ERROR_NULL ( fn, XLAL_EFAILED );
}
/* add duty_X-weighted SXinv to Q */
UINT4 iFreq;
for ( iFreq = 0; iFreq < numFreqs; iFreq ++ )
Q->data->data[iFreq] += duty_X * SXinv->data->data[iFreq] / numSFTsInSeg;
} /* for X < numIFOs */
/* clean up, free memory */
XLALDestroyREAL8FrequencySeries ( SXinv );
return Q;
} /* XLALComputeSegmentDataQ() */
/** Write given REAL8FrequencySeries into file
*/
int
XLALWriteREAL8FrequencySeries_to_file ( const REAL8FrequencySeries *series, /**< [in] frequency-series to write to file */
const char *fname /**< [in] filename to write into */
)
{
const char *fn = __func__;
/* check input consistency */
if ( !series || !fname ) {
XLALPrintError ("%s: invalid NULL input.\n", fn );
XLAL_ERROR ( fn, XLAL_EINVAL );
}
FILE *fp;
if ( ( fp = fopen ( fname, "wb" )) == NULL ) {
XLALPrintError ("%s: failed to open file '%s' for writing.\n", fn, fname );
XLAL_ERROR ( fn, XLAL_ESYS );
}
/* write header info in comments */
if ( XLAL_SUCCESS != XLALOutputVersionString ( fp, 0 ) )
XLAL_ERROR ( fn, XLAL_EFUNC );
fprintf ( fp, "%%%% name = '%s'\n", series->name );
fprintf ( fp, "%%%% epoch = {%d, %d}\n", series->epoch.gpsSeconds, series->epoch.gpsNanoSeconds );
fprintf ( fp, "%%%% f0 = %f Hz\n", series->f0 );
fprintf ( fp, "%%%% deltaF = %g Hz\n", series->deltaF );
CHAR unitStr[1024];
if ( XLALUnitAsString( &unitStr[0], sizeof(unitStr)-1, &series->sampleUnits ) == NULL ) {
XLALPrintError ("%s: XLALUnitAsString() failed with xlalErrno = %d.\n", fn, xlalErrno );
XLAL_ERROR ( fn, XLAL_EFUNC );
}
fprintf ( fp, "%%%% Units = %s\n", unitStr );
fprintf ( fp, "%%%% Freq [Hz] Data(Freq)\n");
UINT4 numBins = series->data->length;
UINT4 iFreq;
for ( iFreq = 0; iFreq < numBins; iFreq ++ )
{
REAL8 thisFreq = series->f0 + iFreq * series->deltaF;
fprintf (fp, "%20.16f %20.16g\n", thisFreq, series->data->data[iFreq] );
}
fclose ( fp );
return XLAL_SUCCESS;
} /* XLALWriteREAL8FrequencySeries_to_file() */
......@@ -929,8 +929,12 @@ int MAIN( int argc, char *argv[]) {
}
REAL8 timeStart = XLALGetTimeOfDay();
/* timing */
REAL8 timeStart = 0.0, timeEnd = 0.0;
REAL8 coherentTime = 0.0, incoherentTime = 0.0, vetoTime = 0.0;
REAL8 timeStamp1 = 0.0, timeStamp2 = 0.0;
if ( uvar_outputTiming )
timeStart = XLALGetTimeOfDay();
/* ################## loop over SKY coarse-grid points ################## */
while(thisScan.state != STATE_FINISHED)
......@@ -1189,6 +1193,10 @@ int MAIN( int argc, char *argv[]) {
u2win = df1dot;
u2winInv = 1.0/u2win; */
/* timing */
if ( uvar_outputTiming )
timeStamp1 = XLALGetTimeOfDay();
/* ----------------------------------------------------------------- */
/************************ Compute F-Statistic ************************/
......@@ -1301,6 +1309,13 @@ int MAIN( int argc, char *argv[]) {
/* --- Holger: This is not needed in U1-only case. Sort the coarse grid in Uindex --- */
/* qsort(coarsegrid.list, (size_t)coarsegrid.length, sizeof(CoarseGridPoint), compareCoarseGridUindex); */
/* timing */
if ( uvar_outputTiming ) {
timeStamp2 = XLALGetTimeOfDay();
coherentTime += timeStamp2 - timeStamp1;
timeStamp1 = timeStamp2;
}
/* ---------- Walk through fine grid and map to coarse grid --------------- */
ifine = 0;
......@@ -1358,6 +1373,13 @@ int MAIN( int argc, char *argv[]) {
fprintf(stderr, " --- Seg: %03d nc_max: %03d avesumTwoFmax: %f \n", k, nc_max, sumTwoFmax/(k+1));
#endif
#endif
/* timing */
if ( uvar_outputTiming ) {
timeStamp2 = XLALGetTimeOfDay();
incoherentTime += timeStamp2 - timeStamp1;
}
} /* end: ------------- MAIN LOOP over Segments --------------------*/
/* ############################################################### */
......@@ -1402,33 +1424,48 @@ int MAIN( int argc, char *argv[]) {
} /* ######## End of while loop over 1st stage SKY coarse-grid points ############ */
/*---------------------------------------------------------------------------------*/
REAL8 timeEnd = XLALGetTimeOfDay();
UINT4 Nrefine = nf1dots_fg;
/* timing */
if ( uvar_outputTiming )
{
FILE *timing_fp = fopen ( uvar_outputTiming, "ab" );
if ( timing_fp == NULL ) {
XLALPrintError ("%s: failed to open timing-file '%s' for appending.\n", __func__, uvar_outputTiming );
return HIERARCHICALSEARCH_EFILE;
}
REAL8 tau = timeEnd - timeStart;
fprintf ( timing_fp, "%d %d %d %d %d %d %d %f\n",
thisScan.numSkyGridPoints, nf1dot, binsFstatSearch, 2 * semiCohPar.extraBinsFstat, nSFTs, nStacks, Nrefine, tau );
fclose ( timing_fp );
}
timeEnd = XLALGetTimeOfDay();
UINT4 Nrefine = nf1dots_fg;
LogPrintf( LOG_NORMAL, "Finished analysis.\n");
/* Also compute F, FX (for line veto statistics) for all candidates in final toplist */
if ( uvar_outputFX ) {
/* timing */
if ( uvar_outputTiming )
timeStamp1 = XLALGetTimeOfDay();
xlalErrno = 0;
XLALComputeExtraStatsForToplist ( semiCohToplist, &stackMultiSFT, &stackMultiNoiseWeights, &stackMultiDetStates, &CFparams, refTimeGPS, tMidGPS, uvar_SignalOnly );
if ( xlalErrno != 0 ) {
XLALPrintError ("%s line %d : XLALComputeLineVetoForToplist() failed with xlalErrno = %d.\n\n", fn, __LINE__, xlalErrno );
return(HIERARCHICALSEARCH_EXLAL);
}
/* timing */
if ( uvar_outputTiming ) {
timeStamp2 = XLALGetTimeOfDay();
vetoTime = timeStamp2 - timeStamp1;
}
}
if ( uvar_outputTiming )
{
FILE *timing_fp = fopen ( uvar_outputTiming, "ab" );
if ( timing_fp == NULL ) {
XLALPrintError ("%s: failed to open timing-file '%s' for appending.\n", __func__, uvar_outputTiming );
return HIERARCHICALSEARCH_EFILE;
}
REAL8 tau = timeEnd - timeStart;
fprintf ( timing_fp, "%d %d %d %d %d %d %d %f %f %f %f\n",
thisScan.numSkyGridPoints, nf1dot, binsFstatSearch, 2 * semiCohPar.extraBinsFstat, nSFTs,
nStacks, Nrefine, tau, coherentTime, incoherentTime, vetoTime );
fclose ( timing_fp );
}
LogPrintf ( LOG_DEBUG, "Writing output ... ");
write_hfs_oputput(uvar_fnameout, semiCohToplist);
LogPrintfVerbatim ( LOG_DEBUG, "done.\n");
......
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