Commit 02f590de authored by Reinhard Prix's avatar Reinhard Prix
Browse files

ComputeFstat and related: pure internal renaming for better clarity

- no API changes
- prepares and makes next API-changing patches simpler
- refs #1936
Original: 51f66fb734cba77950e2601600617abcd76482c0
parent d67c2c5c
...@@ -1361,15 +1361,15 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar ) ...@@ -1361,15 +1361,15 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar )
} /* extrapolate spin-range */ } /* extrapolate spin-range */
/* if single-only flag is given, assume a PSD with sqrt(S) = 1.0 */ /* if single-only flag is given, assume a PSD with sqrt(S) = 1.0 */
MultiNoiseFloor assumeSqrtSX, *p_assumeSqrtSX; MultiNoiseFloor s_assumeSqrtSX, *assumeSqrtSX;
if ( uvar->SignalOnly ) { if ( uvar->SignalOnly ) {
assumeSqrtSX.length = XLALCountIFOsInCatalog(catalog); s_assumeSqrtSX.length = XLALCountIFOsInCatalog(catalog);
for (UINT4 X = 0; X < assumeSqrtSX.length; ++X) { for (UINT4 X = 0; X < s_assumeSqrtSX.length; ++X) {
assumeSqrtSX.sqrtSn[X] = 1.0; s_assumeSqrtSX.sqrtSn[X] = 1.0;
} }
p_assumeSqrtSX = &assumeSqrtSX; assumeSqrtSX = &s_assumeSqrtSX;
} else { } else {
p_assumeSqrtSX = NULL; assumeSqrtSX = NULL;
} }
PulsarParamsVector *injectSources = NULL; PulsarParamsVector *injectSources = NULL;
...@@ -1378,7 +1378,7 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar ) ...@@ -1378,7 +1378,7 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar )
extraParams.Dterms = uvar->Dterms; extraParams.Dterms = uvar->Dterms;
extraParams.SSBprec = uvar->SSBprecision; extraParams.SSBprec = uvar->SSBprecision;
cfg->Fstat_in = XLALCreateFstatInput( catalog, fCoverMin, fCoverMax, cfg->Fstat_in = XLALCreateFstatInput( catalog, fCoverMin, fCoverMax,
injectSources, injectSqrtSX, p_assumeSqrtSX, uvar->RngMedWindow, injectSources, injectSqrtSX, assumeSqrtSX, uvar->RngMedWindow,
cfg->ephemeris, cfg->FstatMethod, &extraParams ); cfg->ephemeris, cfg->FstatMethod, &extraParams );
XLAL_CHECK ( cfg->Fstat_in != NULL, XLAL_EFUNC, "XLALCreateFstatInput() failed with errno=%d", xlalErrno); XLAL_CHECK ( cfg->Fstat_in != NULL, XLAL_EFUNC, "XLALCreateFstatInput() failed with errno=%d", xlalErrno);
XLALDestroySFTCatalog(catalog); XLALDestroySFTCatalog(catalog);
......
...@@ -2038,16 +2038,16 @@ void SetUpSFTs( LALStatus *status, /**< pointer to LALStatus structure */ ...@@ -2038,16 +2038,16 @@ void SetUpSFTs( LALStatus *status, /**< pointer to LALStatus structure */
for (k = 0; k < in->nStacks; k++) { for (k = 0; k < in->nStacks; k++) {
/* if single-only flag is given, assume a PSD with sqrt(S) = 1.0 */ /* if single-only flag is given, assume a PSD with sqrt(S) = 1.0 */
MultiNoiseFloor assumeSqrtSX, *p_assumeSqrtSX; MultiNoiseFloor s_assumeSqrtSX, *assumeSqrtSX;
if ( in->SignalOnly ) { if ( in->SignalOnly ) {
const SFTCatalog *catalog_k = &(catalogSeq.data[k]); const SFTCatalog *catalog_k = &(catalogSeq.data[k]);
assumeSqrtSX.length = XLALCountIFOsInCatalog ( catalog_k ); s_assumeSqrtSX.length = XLALCountIFOsInCatalog ( catalog_k );
for (UINT4 X = 0; X < assumeSqrtSX.length; ++X) { for (UINT4 X = 0; X < s_assumeSqrtSX.length; ++X) {
assumeSqrtSX.sqrtSn[X] = 1.0; s_assumeSqrtSX.sqrtSn[X] = 1.0;
} }
p_assumeSqrtSX = &assumeSqrtSX; assumeSqrtSX = &s_assumeSqrtSX;
} else { } else {
p_assumeSqrtSX = NULL; assumeSqrtSX = NULL;
} }
PulsarParamsVector *injectSources = NULL; PulsarParamsVector *injectSources = NULL;
...@@ -2055,7 +2055,7 @@ void SetUpSFTs( LALStatus *status, /**< pointer to LALStatus structure */ ...@@ -2055,7 +2055,7 @@ void SetUpSFTs( LALStatus *status, /**< pointer to LALStatus structure */
/* ----- create Fstat input data struct ----- */ /* ----- create Fstat input data struct ----- */
(*p_Fstat_in_vec)->data[k] = XLALCreateFstatInput ( &catalogSeq.data[k], freqmin, freqmax, (*p_Fstat_in_vec)->data[k] = XLALCreateFstatInput ( &catalogSeq.data[k], freqmin, freqmax,
injectSources, injectSqrtSX, p_assumeSqrtSX, in->blocksRngMed, injectSources, injectSqrtSX, assumeSqrtSX, in->blocksRngMed,
in->edat, in->Fmethod, &extraParams ); in->edat, in->Fmethod, &extraParams );
if ( (*p_Fstat_in_vec)->data[k] == NULL ) { if ( (*p_Fstat_in_vec)->data[k] == NULL ) {
XLALPrintError("%s: XLALCreateFstatInput() failed with errno=%d", __func__, xlalErrno); XLALPrintError("%s: XLALCreateFstatInput() failed with errno=%d", __func__, xlalErrno);
......
...@@ -42,9 +42,9 @@ ...@@ -42,9 +42,9 @@
// Common input data for F-statistic methods // Common input data for F-statistic methods
typedef struct { typedef struct {
MultiLALDetector detectors; // List of detectors MultiLALDetector detectors; // List of detectors
MultiLIGOTimeGPSVector *timestamps; // Multi-detector list of SFT timestamps MultiLIGOTimeGPSVector *multiTimestamps; // Multi-detector list of SFT timestamps
MultiNoiseWeights *noiseWeights; // Multi-detector noise weights MultiNoiseWeights *multiNoiseWeights; // Multi-detector noise weights
MultiDetectorStateSeries *detectorStates; // Multi-detector state series MultiDetectorStateSeries *multiDetectorStates; // Multi-detector state series
const EphemerisData *ephemerides; // Ephemerides for the time-span of the SFTs const EphemerisData *ephemerides; // Ephemerides for the time-span of the SFTs
SSBprecision SSBprec; // Barycentric transformation precision SSBprecision SSBprec; // Barycentric transformation precision
FstatMethodType FstatMethod; // Method to use for computing the F-statistic FstatMethodType FstatMethod; // Method to use for computing the F-statistic
...@@ -320,7 +320,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT ...@@ -320,7 +320,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
REAL8 minFreqMethod, maxFreqMethod; REAL8 minFreqMethod, maxFreqMethod;
REAL8 minFreqFull, maxFreqFull; REAL8 minFreqFull, maxFreqFull;
{ {
// Determine whether the method being used requires extra frequency bins // Determine whether the method being used requires extra SFT frequency bins
int extraBinsMethod = 0; int extraBinsMethod = 0;
if ( input->demod != NULL ) if ( input->demod != NULL )
{ {
...@@ -359,7 +359,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT ...@@ -359,7 +359,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
// Extract detectors and timestamps from SFTs // Extract detectors and timestamps from SFTs
XLAL_CHECK_NULL ( XLALMultiLALDetectorFromMultiSFTs ( &common->detectors, multiSFTs ) == XLAL_SUCCESS, XLAL_EFUNC ); XLAL_CHECK_NULL ( XLALMultiLALDetectorFromMultiSFTs ( &common->detectors, multiSFTs ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK_NULL ( ( common->timestamps = XLALExtractMultiTimestampsFromSFTs ( multiSFTs ) ) != NULL, XLAL_EFUNC ); XLAL_CHECK_NULL ( ( common->multiTimestamps = XLALExtractMultiTimestampsFromSFTs ( multiSFTs ) ) != NULL, XLAL_EFUNC );
} }
else else
...@@ -370,7 +370,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT ...@@ -370,7 +370,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
// Extract detectors and timestamps from multi-view of SFT catalog // Extract detectors and timestamps from multi-view of SFT catalog
XLAL_CHECK_NULL ( XLALMultiLALDetectorFromMultiSFTCatalogView ( &common->detectors, multiSFTcatalog ) == XLAL_SUCCESS, XLAL_EFUNC ); XLAL_CHECK_NULL ( XLALMultiLALDetectorFromMultiSFTCatalogView ( &common->detectors, multiSFTcatalog ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK_NULL ( ( common->timestamps = XLALTimestampsFromMultiSFTCatalogView ( multiSFTcatalog ) ) != NULL, XLAL_EFUNC ); XLAL_CHECK_NULL ( ( common->multiTimestamps = XLALTimestampsFromMultiSFTCatalogView ( multiSFTcatalog ) ) != NULL, XLAL_EFUNC );
// Cleanup // Cleanup
XLALDestroyMultiSFTCatalogView ( multiSFTcatalog ); XLALDestroyMultiSFTCatalogView ( multiSFTcatalog );
...@@ -388,7 +388,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT ...@@ -388,7 +388,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
MFDparams.fMin = minFreqFull; MFDparams.fMin = minFreqFull;
MFDparams.Band = maxFreqFull - minFreqFull; MFDparams.Band = maxFreqFull - minFreqFull;
MFDparams.multiIFO = common->detectors; MFDparams.multiIFO = common->detectors;
MFDparams.multiTimestamps = *(common->timestamps); MFDparams.multiTimestamps = *(common->multiTimestamps);
MFDparams.randSeed = extraParams->randSeed; MFDparams.randSeed = extraParams->randSeed;
// Set noise floors if sqrtSX is given; otherwise noise floors are zero // Set noise floors if sqrtSX is given; otherwise noise floors are zero
...@@ -422,7 +422,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT ...@@ -422,7 +422,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
XLAL_CHECK_NULL ( (runningMedian = XLALNormalizeMultiSFTVect ( multiSFTs, runningMedianWindow, assumeSqrtSX )) != NULL, XLAL_EFUNC ); XLAL_CHECK_NULL ( (runningMedian = XLALNormalizeMultiSFTVect ( multiSFTs, runningMedianWindow, assumeSqrtSX )) != NULL, XLAL_EFUNC );
// Calculate SFT noise weights from PSD // Calculate SFT noise weights from PSD
XLAL_CHECK_NULL ( (common->noiseWeights = XLALComputeMultiNoiseWeights ( runningMedian, runningMedianWindow, 0 )) != NULL, XLAL_EFUNC ); XLAL_CHECK_NULL ( (common->multiNoiseWeights = XLALComputeMultiNoiseWeights ( runningMedian, runningMedianWindow, 0 )) != NULL, XLAL_EFUNC );
// at this point we're done with running-median noise estimation and can 'trim' the SFTs back to // at this point we're done with running-median noise estimation and can 'trim' the SFTs back to
// the width actually required by the Fstat-methods *methods*. // the width actually required by the Fstat-methods *methods*.
...@@ -432,7 +432,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT ...@@ -432,7 +432,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
// Get detector states, with a timestamp shift of Tsft/2 // Get detector states, with a timestamp shift of Tsft/2
const REAL8 tOffset = 0.5 * Tsft; const REAL8 tOffset = 0.5 * Tsft;
XLAL_CHECK_NULL ( (common->detectorStates = XLALGetMultiDetectorStates ( common->timestamps, &common->detectors, ephemerides, tOffset )) != NULL, XLAL_EFUNC ); XLAL_CHECK_NULL ( (common->multiDetectorStates = XLALGetMultiDetectorStates ( common->multiTimestamps, &common->detectors, ephemerides, tOffset )) != NULL, XLAL_EFUNC );
// Save ephemerides and SSB precision // Save ephemerides and SSB precision
common->ephemerides = ephemerides; common->ephemerides = ephemerides;
...@@ -488,7 +488,7 @@ XLALGetFstatInputTimestamps ( const FstatInput* input ///< [in] \c FstatInput st ...@@ -488,7 +488,7 @@ XLALGetFstatInputTimestamps ( const FstatInput* input ///< [in] \c FstatInput st
XLAL_CHECK_NULL ( input != NULL, XLAL_EINVAL ); XLAL_CHECK_NULL ( input != NULL, XLAL_EINVAL );
XLAL_CHECK_NULL ( input->common != NULL, XLAL_EINVAL, "'input' has not yet been set up" ); XLAL_CHECK_NULL ( input->common != NULL, XLAL_EINVAL, "'input' has not yet been set up" );
return input->common->timestamps; return input->common->multiTimestamps;
} // XLALGetFstatInputTimestamps() } // XLALGetFstatInputTimestamps()
...@@ -503,7 +503,7 @@ XLALGetFstatInputNoiseWeights ( const FstatInput* input ///< [in] \c FstatIn ...@@ -503,7 +503,7 @@ XLALGetFstatInputNoiseWeights ( const FstatInput* input ///< [in] \c FstatIn
XLAL_CHECK_NULL ( input != NULL, XLAL_EINVAL ); XLAL_CHECK_NULL ( input != NULL, XLAL_EINVAL );
XLAL_CHECK_NULL ( input->common != NULL, XLAL_EINVAL, "'input' has not yet been set up" ); XLAL_CHECK_NULL ( input->common != NULL, XLAL_EINVAL, "'input' has not yet been set up" );
return input->common->noiseWeights; return input->common->multiNoiseWeights;
} // XLALGetFstatInputNoiseWeights() } // XLALGetFstatInputNoiseWeights()
...@@ -518,7 +518,7 @@ XLALGetFstatInputDetectorStates ( const FstatInput* input ///< [in] \c FstatInpu ...@@ -518,7 +518,7 @@ XLALGetFstatInputDetectorStates ( const FstatInput* input ///< [in] \c FstatInpu
XLAL_CHECK_NULL ( input != NULL, XLAL_EINVAL ); XLAL_CHECK_NULL ( input != NULL, XLAL_EINVAL );
XLAL_CHECK_NULL ( input->common != NULL, XLAL_EINVAL, "'input' has not yet been set up" ); XLAL_CHECK_NULL ( input->common != NULL, XLAL_EINVAL, "'input' has not yet been set up" );
return input->common->detectorStates; return input->common->multiDetectorStates;
} // XLALGetFstatInputDetectorStates() } // XLALGetFstatInputDetectorStates()
...@@ -667,9 +667,9 @@ XLALDestroyFstatInput ( FstatInput* input ///< [in] \c FstatInput structure to b ...@@ -667,9 +667,9 @@ XLALDestroyFstatInput ( FstatInput* input ///< [in] \c FstatInput structure to b
if (input->common != NULL) if (input->common != NULL)
{ {
XLALDestroyMultiTimestamps ( input->common->timestamps ); XLALDestroyMultiTimestamps ( input->common->multiTimestamps );
XLALDestroyMultiNoiseWeights ( input->common->noiseWeights ); XLALDestroyMultiNoiseWeights ( input->common->multiNoiseWeights );
XLALDestroyMultiDetectorStateSeries ( input->common->detectorStates ); XLALDestroyMultiDetectorStateSeries ( input->common->multiDetectorStates );
XLALFree ( input->common ); XLALFree ( input->common );
} }
if (input->demod != NULL) if (input->demod != NULL)
......
...@@ -117,8 +117,8 @@ ComputeFstat_Demod ( FstatResults* Fstats, ...@@ -117,8 +117,8 @@ ComputeFstat_Demod ( FstatResults* Fstats,
PulsarDopplerParams thisPoint = Fstats->doppler; PulsarDopplerParams thisPoint = Fstats->doppler;
const REAL8 fStart = thisPoint.fkdot[0]; const REAL8 fStart = thisPoint.fkdot[0];
const MultiSFTVector *multiSFTs = demod->multiSFTs; const MultiSFTVector *multiSFTs = demod->multiSFTs;
const MultiNoiseWeights *multiWeights = common->noiseWeights; const MultiNoiseWeights *multiWeights = common->multiNoiseWeights;
const MultiDetectorStateSeries *multiDetStates = common->detectorStates; const MultiDetectorStateSeries *multiDetStates = common->multiDetectorStates;
UINT4 numDetectors = multiSFTs->length; UINT4 numDetectors = multiSFTs->length;
XLAL_CHECK ( multiDetStates->length == numDetectors, XLAL_EINVAL ); XLAL_CHECK ( multiDetStates->length == numDetectors, XLAL_EINVAL );
......
...@@ -208,7 +208,7 @@ SetupFstatInput_Resamp ( FstatInput_Resamp *resamp, ...@@ -208,7 +208,7 @@ SetupFstatInput_Resamp ( FstatInput_Resamp *resamp,
numSamplesInMax = MYMAX ( numSamplesInMax, numSamplesInX ); numSamplesInMax = MYMAX ( numSamplesInMax, numSamplesInX );
XLAL_CHECK ( (resamp->prev_multiTimeSeries_SRC->data[X] = XLALCreateCOMPLEX8TimeSeries ( "", &epoch0, fHet0, dt0, &emptyLALUnit, numSamplesInX )) != NULL, XLAL_EFUNC ); XLAL_CHECK ( (resamp->prev_multiTimeSeries_SRC->data[X] = XLALCreateCOMPLEX8TimeSeries ( "", &epoch0, fHet0, dt0, &emptyLALUnit, numSamplesInX )) != NULL, XLAL_EFUNC );
UINT4 numTimestamps = common->timestamps->data[X]->length; UINT4 numTimestamps = common->multiTimestamps->data[X]->length;
XLAL_CHECK ( (resamp->prev_multiSFTinds_SRC->data[X] = XLALCreateUINT4Vector ( 2 * numTimestamps )) != NULL, XLAL_EFUNC ); XLAL_CHECK ( (resamp->prev_multiSFTinds_SRC->data[X] = XLALCreateUINT4Vector ( 2 * numTimestamps )) != NULL, XLAL_EFUNC );
} // for X < numDetectors } // for X < numDetectors
...@@ -279,7 +279,7 @@ ComputeFstat_Resamp ( FstatResults* Fstats, ...@@ -279,7 +279,7 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
else else
{ {
XLALDestroyMultiAMCoeffs ( resamp->prev_multiAMcoef ); XLALDestroyMultiAMCoeffs ( resamp->prev_multiAMcoef );
XLAL_CHECK ( (multiAMcoef = XLALComputeMultiAMCoeffs ( common->detectorStates, common->noiseWeights, skypos )) != NULL, XLAL_EFUNC ); XLAL_CHECK ( (multiAMcoef = XLALComputeMultiAMCoeffs ( common->multiDetectorStates, common->multiNoiseWeights, skypos )) != NULL, XLAL_EFUNC );
resamp->prev_multiAMcoef = multiAMcoef; resamp->prev_multiAMcoef = multiAMcoef;
} }
...@@ -292,7 +292,7 @@ ComputeFstat_Resamp ( FstatResults* Fstats, ...@@ -292,7 +292,7 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
else else
{ {
XLALDestroyMultiSSBtimes ( resamp->prev_multiTimingSRC ); XLALDestroyMultiSSBtimes ( resamp->prev_multiTimingSRC );
XLAL_CHECK ( (multiTimingSRC = XLALGetMultiSSBtimes ( common->detectorStates, skypos, thisPoint.refTime, common->SSBprec )) != NULL, XLAL_EFUNC ); XLAL_CHECK ( (multiTimingSRC = XLALGetMultiSSBtimes ( common->multiDetectorStates, skypos, thisPoint.refTime, common->SSBprec )) != NULL, XLAL_EFUNC );
if ( thisPoint.asini > 0 ) { if ( thisPoint.asini > 0 ) {
XLAL_CHECK ( XLALAddMultiBinaryTimes ( &multiTimingSRC, multiTimingSRC, &thisPoint ) == XLAL_SUCCESS, XLAL_EFUNC ); XLAL_CHECK ( XLALAddMultiBinaryTimes ( &multiTimingSRC, multiTimingSRC, &thisPoint ) == XLAL_SUCCESS, XLAL_EFUNC );
} }
...@@ -303,7 +303,7 @@ ComputeFstat_Resamp ( FstatResults* Fstats, ...@@ -303,7 +303,7 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
// ----- if NOT same SRC timing OR same frequency-resolution (ie numSamplesOut)? --> recompute SRC-frame timeseries // ----- if NOT same SRC timing OR same frequency-resolution (ie numSamplesOut)? --> recompute SRC-frame timeseries
if ( ! ( same_numSamplesOut && same_skypos && same_refTime && same_binary ) ) if ( ! ( same_numSamplesOut && same_skypos && same_refTime && same_binary ) )
{ {
XLAL_CHECK ( XLALBarycentricResampleMultiCOMPLEX8TimeSeries ( resamp->prev_multiTimeSeries_SRC, resamp->prev_multiSFTinds_SRC, multiTimeSeries_DET, common->timestamps, multiTimingSRC, dt_SRC ) XLAL_CHECK ( XLALBarycentricResampleMultiCOMPLEX8TimeSeries ( resamp->prev_multiTimeSeries_SRC, resamp->prev_multiSFTinds_SRC, multiTimeSeries_DET, common->multiTimestamps, multiTimingSRC, dt_SRC )
== XLAL_SUCCESS, XLAL_EFUNC ); == XLAL_SUCCESS, XLAL_EFUNC );
} }
......
...@@ -89,7 +89,6 @@ main ( int argc, char *argv[] ) ...@@ -89,7 +89,6 @@ main ( int argc, char *argv[] )
fprintf ( stderr, "Tseg = %.1f d, numSegments = %" LAL_INT4_FORMAT ", Freq = %.1f Hz, f1dot = %.1e Hz/s, dFreq = %.1e Hz, numFreqBins = %" LAL_INT4_FORMAT ", FreqBand = %.2f Hz, Tsft = %.0f s\n", fprintf ( stderr, "Tseg = %.1f d, numSegments = %" LAL_INT4_FORMAT ", Freq = %.1f Hz, f1dot = %.1e Hz/s, dFreq = %.1e Hz, numFreqBins = %" LAL_INT4_FORMAT ", FreqBand = %.2f Hz, Tsft = %.0f s\n",
uvar->Tseg / 86400.0, uvar->numSegments, uvar->Freq, uvar->f1dot, dFreq, uvar->numFreqBins, FreqBand, uvar->Tsft ); uvar->Tseg / 86400.0, uvar->numSegments, uvar->Freq, uvar->f1dot, dFreq, uvar->numFreqBins, FreqBand, uvar->Tsft );
// ---------- end: handle user input ---------- // ---------- end: handle user input ----------
EphemerisData *ephem; EphemerisData *ephem;
REAL8 maxMem0 = XLALGetPeakHeapUsageMB(); REAL8 maxMem0 = XLALGetPeakHeapUsageMB();
XLAL_CHECK ( (ephem = XLALInitBarycenter ( TEST_DATA_DIR "earth00-19-DE405.dat.gz", TEST_DATA_DIR "sun00-19-DE405.dat.gz" )) != NULL, XLAL_EFUNC ); XLAL_CHECK ( (ephem = XLALInitBarycenter ( TEST_DATA_DIR "earth00-19-DE405.dat.gz", TEST_DATA_DIR "sun00-19-DE405.dat.gz" )) != NULL, XLAL_EFUNC );
...@@ -99,18 +98,9 @@ main ( int argc, char *argv[] ) ...@@ -99,18 +98,9 @@ main ( int argc, char *argv[] )
XLALPrintInfo ("mem(ephemeris) = %.1f MB\n", memEphem ); XLALPrintInfo ("mem(ephemeris) = %.1f MB\n", memEphem );
// ----- setup injection and data parameters // ----- setup injection and data parameters
UINT4 numDetectors = 2;
// use signal-only injections to avoid usage of different noise bins to contaminate error-comparison
MultiNoiseFloor XLAL_INIT_DECL(injectNoiseFloor);
injectNoiseFloor.length = numDetectors;
injectNoiseFloor.sqrtSn[0] = 1;
injectNoiseFloor.sqrtSn[1] = 2;
LALStringVector *detNames = NULL; LALStringVector *detNames = NULL;
XLAL_CHECK ( (detNames = XLALCreateStringVector ( "H1", "L1", NULL )) != NULL, XLAL_EFUNC ); XLAL_CHECK ( (detNames = XLALCreateStringVector ( "H1", "L1", NULL )) != NULL, XLAL_EFUNC );
MultiLALDetector XLAL_INIT_DECL(detInfo); UINT4 numDetectors = detNames->length;
XLAL_CHECK ( XLALParseMultiLALDetector ( &detInfo, detNames ) == XLAL_SUCCESS, XLAL_EFUNC );
LIGOTimeGPS startTime = {711595934, 0}; LIGOTimeGPS startTime = {711595934, 0};
LIGOTimeGPS startTime_l = startTime; LIGOTimeGPS startTime_l = startTime;
...@@ -163,16 +153,21 @@ main ( int argc, char *argv[] ) ...@@ -163,16 +153,21 @@ main ( int argc, char *argv[] )
extraParams.SSBprec = SSBPREC_RELATIVISTICOPT; extraParams.SSBprec = SSBPREC_RELATIVISTICOPT;
extraParams.Dterms = Dterms; // constant value that works for all Demod methods extraParams.Dterms = Dterms; // constant value that works for all Demod methods
// ----- prepare input data with injection for all available methods MultiNoiseFloor XLAL_INIT_DECL(injectSqrtSX);
FstatQuantities whatToCompute = (FSTATQ_2F | FSTATQ_2F_PER_DET); injectSqrtSX.length = numDetectors;
for ( UINT4 X=0; X < numDetectors; X ++ ) {
injectSqrtSX.sqrtSn[X] = 1;
}
FstatInput **inputs; FstatInputVector *inputs;
XLAL_CHECK ( (inputs = XLALCalloc ( uvar->numSegments, sizeof(inputs[0]))) != NULL, XLAL_ENOMEM ); XLAL_CHECK ( (inputs = XLALCreateFstatInputVector ( uvar->numSegments )) != NULL, XLAL_EFUNC );
for ( INT4 l = 0; l < uvar->numSegments; l ++ ) for ( INT4 l = 0; l < uvar->numSegments; l ++ )
{ {
XLAL_CHECK ( (inputs[l] = XLALCreateFstatInput ( catalogs[l], minCoverFreq, maxCoverFreq, NULL, &injectNoiseFloor, NULL, rngMed, ephem, FstatMethod, &extraParams )) != NULL, XLAL_EFUNC ); XLAL_CHECK ( (inputs->data[l] = XLALCreateFstatInput ( catalogs[l], minCoverFreq, maxCoverFreq, NULL, &injectSqrtSX, NULL, rngMed, ephem, FstatMethod, &extraParams )) != NULL, XLAL_EFUNC );
} }
// REAL8 memMaxSetup = XLALGetPeakHeapUsageMB() - memBase;
// ----- compute Fstatistics over segments
FstatQuantities whatToCompute = (FSTATQ_2F | FSTATQ_2F_PER_DET);
FstatResults **results; FstatResults **results;
XLAL_CHECK ( (results = XLALCalloc ( uvar->numSegments, sizeof(results[0]))) != NULL, XLAL_ENOMEM ); XLAL_CHECK ( (results = XLALCalloc ( uvar->numSegments, sizeof(results[0]))) != NULL, XLAL_ENOMEM );
...@@ -181,12 +176,12 @@ main ( int argc, char *argv[] ) ...@@ -181,12 +176,12 @@ main ( int argc, char *argv[] )
{ {
// call it once to initialize buffering, don't count this time // call it once to initialize buffering, don't count this time
REAL8 tic = XLALGetTimeOfDay(); REAL8 tic = XLALGetTimeOfDay();
XLAL_CHECK ( XLALComputeFstat ( &results[l], inputs[l], &Doppler, dFreq, uvar->numFreqBins, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC ); XLAL_CHECK ( XLALComputeFstat ( &results[l], inputs->data[l], &Doppler, dFreq, uvar->numFreqBins, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
REAL8 toc = XLALGetTimeOfDay(); REAL8 toc = XLALGetTimeOfDay();
tauFSumUnbuffered += ( toc - tic ); tauFSumUnbuffered += ( toc - tic );
// now call it with full buffering to get converged runtime per template (assuming many templates per skypoint or per binary params) // now call it with full buffering to get converged runtime per template (assuming many templates per skypoint or per binary params)
tic = XLALGetTimeOfDay(); tic = XLALGetTimeOfDay();
XLAL_CHECK ( XLALComputeFstat ( &results[l], inputs[l], &Doppler, dFreq, uvar->numFreqBins, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC ); XLAL_CHECK ( XLALComputeFstat ( &results[l], inputs->data[l], &Doppler, dFreq, uvar->numFreqBins, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
toc = XLALGetTimeOfDay(); toc = XLALGetTimeOfDay();
tauFSumBuffered += (toc - tic); tauFSumBuffered += (toc - tic);
} // for l < numSegments } // for l < numSegments
...@@ -197,14 +192,14 @@ main ( int argc, char *argv[] ) ...@@ -197,14 +192,14 @@ main ( int argc, char *argv[] )
fprintf (stderr, "%-15s: tauF1Buffered = %.1e s (=>tauF0SFT = %1.e s), tauF1Unbuffered = %.1e s, memSFTs = %.1f MB, memMaxCompute = %.1f MB\n", fprintf (stderr, "%-15s: tauF1Buffered = %.1e s (=>tauF0SFT = %1.e s), tauF1Unbuffered = %.1e s, memSFTs = %.1f MB, memMaxCompute = %.1f MB\n",
XLALGetFstatMethodName ( FstatMethod ), tauF1Buffered, numDetectors * tauF1Buffered / numSFTsPerSeg, tauF1Unbuffered, memSFTs, memMaxCompute ); XLALGetFstatMethodName ( FstatMethod ), tauF1Buffered, numDetectors * tauF1Buffered / numSFTsPerSeg, tauF1Unbuffered, memSFTs, memMaxCompute );
// ----- free memory ----------
XLALDestroyFstatInputVector ( inputs );
for ( INT4 l = 0; l < uvar->numSegments; l ++ ) for ( INT4 l = 0; l < uvar->numSegments; l ++ )
{ {
XLALDestroySFTCatalog ( catalogs[l] ); XLALDestroySFTCatalog ( catalogs[l] );
XLALDestroyFstatInput ( inputs[l] );
XLALDestroyFstatResults ( results[l] ); XLALDestroyFstatResults ( results[l] );
} }
XLALFree ( catalogs ); XLALFree ( catalogs );
XLALFree ( inputs );
XLALFree ( results ); XLALFree ( results );
XLALDestroyUserVars(); XLALDestroyUserVars();
......
...@@ -42,37 +42,33 @@ main ( int argc, char *argv[] ) ...@@ -42,37 +42,33 @@ main ( int argc, char *argv[] )
XLAL_CHECK ( argc == 1, XLAL_EINVAL, "No input arguments allowed.\n" ); XLAL_CHECK ( argc == 1, XLAL_EINVAL, "No input arguments allowed.\n" );
XLAL_CHECK ( argv != NULL, XLAL_EINVAL ); XLAL_CHECK ( argv != NULL, XLAL_EINVAL );
// ----- load ephemeris
EphemerisData *ephem; EphemerisData *ephem;
XLAL_CHECK ( (ephem = XLALInitBarycenter ( TEST_DATA_DIR "earth00-19-DE405.dat.gz", TEST_DATA_DIR "sun00-19-DE405.dat.gz" )) != NULL, XLAL_EFUNC ); XLAL_CHECK ( (ephem = XLALInitBarycenter ( TEST_DATA_DIR "earth00-19-DE405.dat.gz", TEST_DATA_DIR "sun00-19-DE405.dat.gz" )) != NULL, XLAL_EFUNC );
// ----- setup injection and data parameters // ----- setup injection and data parameters
UINT4 numDetectors = 2;
// use signal-only injections to avoid usage of different noise bins to contaminate error-comparison
MultiNoiseFloor XLAL_INIT_DECL(injectNoiseFloor);
injectNoiseFloor.length = numDetectors;
injectNoiseFloor.sqrtSn[0] = 1;
injectNoiseFloor.sqrtSn[1] = 3;
MultiNoiseFloor XLAL_INIT_DECL(assumeNoiseFloor);
assumeNoiseFloor.length = numDetectors;
assumeNoiseFloor.sqrtSn[0] = 1;
assumeNoiseFloor.sqrtSn[1] = 3;
LALStringVector *detNames = NULL; LALStringVector *detNames = NULL;
XLAL_CHECK ( (detNames = XLALCreateStringVector ( "H1", "L1", NULL )) != NULL, XLAL_EFUNC ); XLAL_CHECK ( (detNames = XLALCreateStringVector ( "H1", "L1", NULL )) != NULL, XLAL_EFUNC );
MultiLALDetector XLAL_INIT_DECL(detInfo); UINT4 numDetectors = detNames->length;
XLAL_CHECK ( XLALParseMultiLALDetector ( &detInfo, detNames ) == XLAL_SUCCESS, XLAL_EFUNC );
// generate and assume some gaussian noise floors
MultiNoiseFloor XLAL_INIT_DECL(injectSqrtSX);
MultiNoiseFloor XLAL_INIT_DECL(assumeSqrtSX);
injectSqrtSX.length = numDetectors;
assumeSqrtSX.length = numDetectors;
for ( UINT4 X = 0; X < numDetectors; X ++ )
{
injectSqrtSX.sqrtSn[X] = 1.0 + 2.0*X;
assumeSqrtSX.sqrtSn[X] = 1.0 + 2.0*X;
}
LIGOTimeGPS startTime = {711595934, 0}; LIGOTimeGPS startTime = {711595934, 0};
REAL8 Tspan = 200 * 3600; REAL8 Tspan = 20 * 3600;
LIGOTimeGPS endTime = startTime; LIGOTimeGPS endTime = startTime;
XLALGPSAdd( &endTime, Tspan ); XLALGPSAdd( &endTime, Tspan );
REAL8 Tsft = 1800; REAL8 Tsft = 1800;
LIGOTimeGPS refTime = { startTime.gpsSeconds - 2.3 * Tspan, 0 }; // reftime in middle of segment LIGOTimeGPS refTime = { startTime.gpsSeconds - 2.3 * Tspan, 0 };
MultiLIGOTimeGPSVector *multiTimestamps; MultiLIGOTimeGPSVector *multiTimestamps;
XLAL_CHECK ( ( multiTimestamps = XLALCalloc ( 1, sizeof(*multiTimestamps))) != NULL, XLAL_ENOMEM ); XLAL_CHECK ( ( multiTimestamps = XLALCalloc ( 1, sizeof(*multiTimestamps))) != NULL, XLAL_ENOMEM );
...@@ -99,13 +95,13 @@ main ( int argc, char *argv[] ) ...@@ -99,13 +95,13 @@ main ( int argc, char *argv[] )
// ----- CW sources to injet ---------- // ----- CW sources to injet ----------
REAL8 Freq = 100.0; REAL8 Freq = 100.0;
PulsarParamsVector *sources; PulsarParamsVector *injectSources;
XLAL_CHECK ( (sources = XLALCreatePulsarParamsVector(1)) != NULL, XLAL_EFUNC ); XLAL_CHECK ( (injectSources = XLALCreatePulsarParamsVector(1)) != NULL, XLAL_EFUNC );
sources->data[0].Amp.h0 = 1; injectSources->data[0].Amp.h0 = 1;
sources->data[0].Amp.cosi = 0.5; injectSources->data[0].Amp.cosi = 0.5;
sources->data[0].Amp.psi = 0.1; injectSources->data[0].Amp.psi = 0.1;
sources->data[0].Amp.phi0 = 1.2; injectSources->data[0].Amp.phi0 = 1.2;
REAL8 asini = 0; // 1.4; // sco-X1 like REAL8 asini = 0; // 1.4; // sco-X1 like
REAL8 Period = 0; // 19 * 3600;// sco-X1 like REAL8 Period = 0; // 19 * 3600;// sco-X1 like
...@@ -123,7 +119,7 @@ main ( int argc, char *argv[] ) ...@@ -123,7 +119,7 @@ main ( int argc, char *argv[] )
Doppler.period = Period; Doppler.period = Period;
Doppler.argp = 0.5; Doppler.argp = 0.5;
sources->data[0].Doppler = Doppler; injectSources->data[0].Doppler = Doppler;
REAL8 dFreq = 0.1 / Tspan; // 10x finer than native FFT resolution REAL8 dFreq = 0.1 / Tspan; // 10x finer than native FFT resolution
REAL8 mis = 0.5; REAL8 mis = 0.5;
...@@ -139,7 +135,7 @@ main ( int argc, char *argv[] ) ...@@ -139,7 +135,7 @@ main ( int argc, char *argv[] )
PulsarSpinRange XLAL_INIT_DECL(spinRange); PulsarSpinRange XLAL_INIT_DECL(spinRange);
spinRange.refTime = refTime; spinRange.refTime = refTime;
memcpy ( &spinRange.fkdot, &sources->data[0].Doppler.fkdot, sizeof(spinRange.fkdot) ); memcpy ( &spinRange.fkdot, &injectSources->data[0].Doppler.fkdot, sizeof(spinRange.fkdot) );
spinRange.fkdotBand[0] = (numFreqBins - 1)*dFreq - 10*LAL_REAL8_EPS; spinRange.fkdotBand[0] = (numFreqBins - 1)*dFreq - 10*LAL_REAL8_EPS;
spinRange.fkdotBand[1] = (numf1dotPoints - 1)*df1dot - 10*LAL_REAL8_EPS; spinRange.fkdotBand[1] = (numf1dotPoints - 1)*df1dot - 10*LAL_REAL8_EPS;
...@@ -163,7 +159,7 @@ main ( int argc, char *argv[] ) ...@@ -163,7 +159,7 @@ main ( int argc, char *argv[] )
if ( !XLALFstatMethodIsAvailable(iMethod) ) { if ( !XLALFstatMethodIsAvailable(iMethod) ) {
continue; continue;
} }
input[iMethod] = XLALCreateFstatInput ( catalog, minCoverFreq, maxCoverFreq, sources, &injectNoiseFloor, &assumeNoiseFloor, 50, ephem, iMethod, &extraParams ); input[iMethod] = XLALCreateFstatInput ( catalog, minCoverFreq, maxCoverFreq, injectSources, &injectSqrtSX, &assumeSqrtSX, 50, ephem, iMethod, &extraParams );
XLAL_CHECK ( input[iMethod] != NULL, XLAL_EFUNC ); XLAL_CHECK ( input[iMethod] != NULL, XLAL_EFUNC );
} }
...@@ -244,11 +240,12 @@ main ( int argc, char *argv[] ) ...@@ -244,11 +240,12 @@ main ( int argc, char *argv[] )
XLALDestroyFstatResults ( results[iMethod] ); XLALDestroyFstatResults ( results[iMethod] );
} // for i < FMETHOD_END } // for i < FMETHOD_END
XLALDestroyPulsarParamsVector ( sources ); XLALDestroyPulsarParamsVector ( injectSources );
XLALDestroySFTCatalog ( catalog ); XLALDestroySFTCatalog ( catalog );
XLALDestroyMultiTimestamps ( multiTimestamps ); XLALDestroyMultiTimestamps ( multiTimestamps );
XLALDestroyStringVector ( detNames ); XLALDestroyStringVector ( detNames );
XLALDestroyEphemerisData ( ephem ); XLALDestroyEphemerisData ( ephem );
LALCheckMemoryLeaks(); LALCheckMemoryLeaks();