Commit 73163496 authored by Reinhard Prix's avatar Reinhard Prix

ComputeFstat_Resamp: more detailed and fine-grained stats and timing info

- sneak this info out via the resampling 'workspace'
- added a new function to append this in readily-readable + plottable form into a file
- refer all timings to single-IFO case
- start using this in CFSBenchmark code
- refs #1954
Original: e98288919993d0b1ff946e9e14c3a7c52a529649
parent 3b0aa642
......@@ -320,7 +320,6 @@ SWIGLAL(INOUT_STRUCTS(FstatResults**, Fstats));
int XLALComputeFstat ( FstatResults **Fstats, FstatInput *input, const PulsarDopplerParams *doppler,
const UINT4 numFreqBins, const FstatQuantities whatToCompute );
FstatWorkspace *XLALGetSharedFstatWorkspace ( FstatInput *input );
void XLALDestroyFstatWorkspace ( FstatWorkspace *ws );
......
......@@ -33,6 +33,7 @@
#define MYMIN(x,y) ( (x) < (y) ? (x) : (y) )
// ----- local constants
#define COLLECT_TIMING 1
#define NUM_FACT 7
static const REAL8 inv_fact[NUM_FACT] = { 1.0, 1.0, (1.0/2.0), (1.0/6.0), (1.0/24.0), (1.0/120.0), (1.0/720.0) };
......@@ -44,6 +45,21 @@ typedef struct tagMultiUINT4Vector
} MultiUINT4Vector;
// ----- workspace ----------
typedef struct tagResampTimingInfo
{ // NOTE: all times refer to a single-detector timing case
REAL8 tauTotal; // total time spent in ComputeFstat_Resamp()
REAL8 tauBary; // time spent in barycentric resampling
REAL8 tauSpin; // time spent in spindown+frequency correction
REAL8 tauAM; // time spent in applying AM coefficients
REAL8 tauFFT; // time spent in FFT
REAL8 tauNorm; // time spent normalizing the final Fa,Fb
REAL8 tauFab2F; // time to compute Fstat from {Fa,Fb}
REAL8 tauMem; // time to realloc and memset-0 arrays
REAL8 tauSumFabX; // time to sum_X Fab^X
REAL8 tauF1Buf; // Resampling timing 'constant': Fstat time per template per detector for a 'buffered' case (same skypos, same numFreqBins)
REAL8 tauF1NoBuf; // Resampling timing 'constant': Fstat time per template per detector for an 'unbuffered' usage (different skypos and numFreqBins)
} ResampTimingInfo;
struct tagFstatWorkspace
{
COMPLEX8Vector *TimeSeriesSpinCorr_SRC; // single-detector SRC-frame spindown-corrected timeseries
......@@ -59,6 +75,8 @@ struct tagFstatWorkspace
COMPLEX8 *FbX_k; // properly normalized F_b^X(f_k) over output bins
COMPLEX8 *Fa_k; // properly normalized F_a(f_k) over output bins
COMPLEX8 *Fb_k; // properly normalized F_b(f_k) over output bins
ResampTimingInfo timingInfo; // temporary storage for collecting timing data
};
struct tagFstatInput_Resamp
......@@ -119,6 +137,8 @@ XLALComputeFaFb_Resamp ( FstatWorkspace *ws,
const AMCoeffs *ab
);
int XLALAppendResampInfo2File ( FILE *fp, const FstatInput *input );
// ==================== function definitions ====================
// ---------- exported API functions ----------
......@@ -169,6 +189,39 @@ XLALDestroyFstatWorkspace ( FstatWorkspace *ws )
} // XLALDestroyFstatWorkspace()
/// debug/optimizer helper function: dump internal info from resampling code into a file
/// if called with input==NULL, output a header-comment-line
int
XLALAppendResampInfo2File ( FILE *fp, const FstatInput *input )
{
XLAL_CHECK ( fp != NULL, XLAL_EINVAL );
if ( input == NULL ) {
fprintf (fp, "%%%%%8s %10s %6s %10s %10s ",
"Nfreq", "NsFFT", "Nsft0", "Ns_DET0", "Ns_SRC0" );
fprintf (fp, "%10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s\n",
"tauTotal", "tauFFT", "tauBary", "tauSpin", "tauAM", "tauNorm", "tauFab2F", "tauMem", "tauSumFabX", "tauF1NoBuf", "tauF1Buf" );
return XLAL_SUCCESS;
}
XLAL_CHECK ( input->resamp != NULL, XLAL_EINVAL );
const FstatInput_Resamp *resamp = input->resamp;
const FstatWorkspace *ws = resamp->ws;
fprintf (fp, "%10d %10d", ws->numFreqBinsOut, ws->numSamplesFFT );
UINT4 numSamples_DETX0 = resamp->multiTimeSeries_DET->data[0]->data->length;
UINT4 numSFTs_X0 = (resamp->prev_multiAMcoef != NULL) ? resamp->prev_multiAMcoef->data[0]->a->length : 0;
COMPLEX8TimeSeries *ts_SRCX0 = resamp->prev_multiTimeSeries_SRC->data[0];
UINT4 numSamples_SRCX0 = ts_SRCX0->data->length;
fprintf (fp, " %6d %10d %10d ", numSFTs_X0, numSamples_DETX0, numSamples_SRCX0 );
const ResampTimingInfo *ti = &(ws->timingInfo);
fprintf (fp, "%10.1e %10.1e %10.1e %10.1e %10.1e %10.1e %10.1e %10.1e %10.1e %10.1e %10.1e\n",
ti->tauTotal, ti->tauFFT, ti->tauBary, ti->tauSpin, ti->tauAM, ti->tauNorm, ti->tauFab2F, ti->tauMem, ti->tauSumFabX, ti->tauF1NoBuf, ti->tauF1Buf );
return XLAL_SUCCESS;
} // XLALAppendResampInfo2File()
// ---------- internal functions ----------
static void
XLALDestroyMultiUINT4Vector ( MultiUINT4Vector *v)
......@@ -312,6 +365,15 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
const FstatQuantities whatToCompute = Fstats->whatWasComputed;
XLAL_CHECK ( !(whatToCompute & FSTATQ_ATOMS_PER_DET), XLAL_EINVAL, "Resampling does not currently support atoms per detector" );
#ifdef COLLECT_TIMING
// collect internal timing info
XLAL_INIT_MEM ( resamp->ws->timingInfo );
ResampTimingInfo *ti = &(resamp->ws->timingInfo);
REAL8 ticStart,tocEnd;
ticStart = XLALGetCPUTime();
REAL8 tic,toc;
#endif
// ----- handy shortcuts ----------
PulsarDopplerParams thisPoint = Fstats->doppler;
const MultiCOMPLEX8TimeSeries *multiTimeSeries_DET = resamp->multiTimeSeries_DET;
......@@ -333,6 +395,10 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
skypos.latitude = thisPoint.Delta;
FstatWorkspace *ws = resamp->ws;
#ifdef COLLECT_TIMING
tic = XLALGetCPUTime();
#endif
// ----- same skyposition? --> reuse antenna-patterns
MultiAMCoeffs *multiAMcoef;
if ( same_skypos && (resamp->prev_multiAMcoef != NULL) )
......@@ -358,6 +424,11 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
== XLAL_SUCCESS, XLAL_EFUNC );
XLALDestroyMultiSSBtimes ( multiTimingSRC );
}
#ifdef COLLECT_TIMING
toc = XLALGetCPUTime();
ti->tauBary = (toc-tic);
#endif
resamp->prev_doppler = thisPoint;
MultiCOMPLEX8TimeSeries *multiTimeSeries_SRC = resamp->prev_multiTimeSeries_SRC;
......@@ -369,6 +440,10 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
UINT4 numFreqBins = Fstats->numFreqBins;
ws->numFreqBinsOut = numFreqBins;
#ifdef COLLECT_TIMING
tic = XLALGetCPUTime();
#endif
// NOTE: we try to use as much existing memory as possible in FstatResults, so we only
// allocate local 'workspace' storage in case there's not already a vector allocated in FstatResults for it
// this also avoid having to copy these results in case the user asked for them to be returned
......@@ -411,6 +486,11 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
memset ( ws->Fa_k, 0, numFreqBins * sizeof(COMPLEX8) );
memset ( ws->Fb_k, 0, numFreqBins * sizeof(COMPLEX8) );
#ifdef COLLECT_TIMING
toc = XLALGetCPUTime();
ti->tauMem = (toc-tic); // this one doesn't scale with number of detector!
#endif
// loop over detectors
for ( UINT4 X=0; X < numDetectors; X++ )
{
......@@ -427,12 +507,19 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
// compute {Fa^X(f_k), Fb^X(f_k)}: results returned via workspace resamp->ws
XLAL_CHECK ( XLALComputeFaFb_Resamp ( resamp->ws, thisPoint, common->dFreq, TimeSeriesX_SRC, SFTindsX_SRC, abX ) == XLAL_SUCCESS, XLAL_EFUNC );
#ifdef COLLECT_TIMING
tic = XLALGetCPUTime();
#endif
for ( UINT4 k = 0; k < numFreqBins; k++ )
{
ws->Fa_k[k] += ws->FaX_k[k];
ws->Fb_k[k] += ws->FbX_k[k];
}
#ifdef COLLECT_TIMING
toc = XLALGetCPUTime();
ti->tauSumFabX += (toc-tic);
tic = toc;
#endif
// ----- if requested: compute per-detector Fstat_X_k
if ( whatToCompute & FSTATQ_2F_PER_DET )
{
......@@ -445,11 +532,20 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
for ( UINT4 k = 0; k < numFreqBins; k ++ )
{
Fstats->twoFPerDet[X][k] = XLALComputeFstatFromFaFb ( ws->FaX_k[k], ws->FbX_k[k], AdX, BdX, CdX, EdX, DdX_inv );
}
} // for k < numFreqBins
} // for k < numFreqBins
} // end: if compute F_X
#ifdef COLLECT_TIMING
toc = XLALGetCPUTime();
ti->tauFab2F += ( toc - tic );
#endif
} // for X < numDetectors
#ifdef COLLECT_TIMING
ti->tauSumFabX /= numDetectors;
ti->tauFab2F /= numDetectors;
tic = XLALGetCPUTime();
#endif
if ( whatToCompute & FSTATQ_2F )
{
for ( UINT4 k=0; k < numFreqBins; k++ )
......@@ -457,6 +553,10 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
Fstats->twoF[k] = XLALComputeFstatFromFaFb ( ws->Fa_k[k], ws->Fb_k[k], Ad, Bd, Cd, Ed, Dd_inv );
} // for k < numFreqBins
} // if FSTATQ_2F
#ifdef COLLECT_TIMING
toc = XLALGetCPUTime();
ti->tauFab2F += ( toc - tic );
#endif
// Return F-atoms per detector
if (whatToCompute & FSTATQ_ATOMS_PER_DET) {
......@@ -481,6 +581,23 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
ws->FbX_k = NULL;
}
#ifdef COLLECT_TIMING
// timings are per-detector
tocEnd = XLALGetCPUTime();
ti->tauTotal = (tocEnd - ticStart);
// rescale all relevant timings to single-IFO case
ti->tauTotal /= numDetectors;
ti->tauBary /= numDetectors;
ti->tauSpin /= numDetectors;
ti->tauAM /= numDetectors;
ti->tauFFT /= numDetectors;
ti->tauNorm /= numDetectors;
// compute 'fundamental' timing numbers per template per detector
ti->tauF1NoBuf = ti->tauTotal / numFreqBins;
ti->tauF1Buf = (ti->tauTotal - ti->tauBary - ti->tauMem) / numFreqBins;
#endif
return XLAL_SUCCESS;
} // ComputeFstat_Resamp()
......@@ -523,15 +640,33 @@ XLALComputeFaFb_Resamp ( FstatWorkspace *ws, //!< [in,out] pre-allocated 'wor
REAL8 fMinFFT = fHet + freqShift - dFreq * NnegBins;
UINT4 offset_bins = (UINT4) lround ( ( FreqOut0 - fMinFFT ) / dFreq );
#ifdef COLLECT_TIMING
// collect some internal timing info
ResampTimingInfo *ti = &(ws->timingInfo);
REAL8 tic,toc;
tic = XLALGetCPUTime();
#endif
// apply spindown phase-factors, store result in 'workspace'
memset ( ws->TimeSeriesSpinCorr_SRC->data, 0, ws->TimeSeriesSpinCorr_SRC->length * sizeof(COMPLEX8));
XLAL_CHECK ( XLALApplySpindownAndFreqShift ( ws->TimeSeriesSpinCorr_SRC, TimeSeries_SRC, SFTinds_SRC, &thisPoint, freqShift ) == XLAL_SUCCESS, XLAL_EFUNC );
#ifdef COLLECT_TIMING
toc = XLALGetCPUTime();
ti->tauSpin += ( toc - tic);
tic = toc;
#endif
// ----- compute FaX_k
// apply amplitude modulation factors {a,b}, store result in zero-padded timeseries for FFTing
memset ( ws->FabX_Raw, 0, ws->numSamplesFFT * sizeof(ws->FabX_Raw[0]) );
XLAL_CHECK ( XLALApplyAmplitudeModulation ( ws->FabX_Raw, ws->TimeSeriesSpinCorr_SRC, SFTinds_SRC, ab->a ) == XLAL_SUCCESS, XLAL_EFUNC );
#ifdef COLLECT_TIMING
toc = XLALGetCPUTime();
ti->tauAM += ( toc - tic);
tic = toc;
#endif
// Fourier transform the resampled Fa(t)
fftwf_execute ( ws->fftplan );
......@@ -544,11 +679,23 @@ XLALComputeFaFb_Resamp ( FstatWorkspace *ws, //!< [in,out] pre-allocated 'wor
} // for k < numFreqBinsOut
#ifdef COLLECT_TIMING
toc = XLALGetCPUTime();
ti->tauFFT += ( toc - tic);
tic = toc;
#endif
// ----- compute FbX_k
// apply amplitude modulation factors {a,b}, store result in zero-padded timeseries for FFTing
memset ( ws->FabX_Raw, 0, ws->numSamplesFFT * sizeof(ws->FabX_Raw[0]) );
XLAL_CHECK ( XLALApplyAmplitudeModulation ( ws->FabX_Raw, ws->TimeSeriesSpinCorr_SRC, SFTinds_SRC, ab->b ) == XLAL_SUCCESS, XLAL_EFUNC );
#ifdef COLLECT_TIMING
toc = XLALGetCPUTime();
ti->tauAM += ( toc - tic);
tic = toc;
#endif
// Fourier transform the resampled Fa(t)
fftwf_execute ( ws->fftplan );
......@@ -561,6 +708,11 @@ XLALComputeFaFb_Resamp ( FstatWorkspace *ws, //!< [in,out] pre-allocated 'wor
} // for k < numFreqBinsOut
#ifdef COLLECT_TIMING
toc = XLALGetCPUTime();
ti->tauFFT += ( toc - tic);
tic = toc;
#endif
// ----- normalization factors to be applied to Fa and Fb:
const REAL8 dtauX = XLALGPSDiff ( &TimeSeries_SRC->epoch, &thisPoint.refTime );
......@@ -575,6 +727,11 @@ XLALComputeFaFb_Resamp ( FstatWorkspace *ws, //!< [in,out] pre-allocated 'wor
ws->FbX_k[k] *= normX_k;
} // for k < numFreqBinsOut
#ifdef COLLECT_TIMING
toc = XLALGetCPUTime();
ti->tauNorm += ( toc - tic);
tic = toc;
#endif
return XLAL_SUCCESS;
......
......@@ -43,8 +43,11 @@ typedef struct
REAL8 Tseg;
REAL8 Tsft;
INT4 numSegments;
CHAR *outputInfo;
} UserInput_t;
int XLALAppendResampInfo2File ( FILE *fp, const FstatInput *input );
// ---------- main ----------
int
main ( int argc, char *argv[] )
......@@ -61,16 +64,22 @@ main ( int argc, char *argv[] )
uvar->Tseg = 60 * 3600;
uvar->numSegments = 90;
uvar->Tsft = 1800;
uvar->outputInfo = XLALStringDuplicate ( "ComputeFstatBenchmark.info" );
XLAL_CHECK ( XLALRegisterUvarMember ( help, BOOLEAN, 'h', HELP, "Print help message" ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALRegisterUvarMember ( FstatMethod, STRING, 0, OPTIONAL, XLALFstatMethodHelpString() ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALRegisterUvarMember ( Freq, REAL8, 0, OPTIONAL, "Search frequency in Hz" ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALRegisterUvarMember ( f1dot, REAL8, 0, OPTIONAL, "Search spindown f1dot in Hz/s" ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALRegisterUvarMember ( FreqResolution, REAL8, 0, OPTIONAL, "Frequency resolution factor 'r' such that dFreq = 1/(r*T)" ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALRegisterUvarMember ( Tseg, REAL8, 0, OPTIONAL, "Coherent segment length" ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALRegisterUvarMember ( numSegments, INT4, 0, OPTIONAL, "Number of semi-coherent segment" ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALRegisterUvarMember ( numFreqBins, INT4Vector, 0, OPTIONAL, "Range of number of frequency bins to search [2-number range input]" ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALRegisterUvarMember ( IFOs, STRINGVector, 0, OPTIONAL, "IFOs to use" ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALRegisterUvarMember ( numTrials, INT4, 0, OPTIONAL, "Number of repeated trials to run (with potentially randomized parameters)" ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALRegisterUvarMember ( outputInfo, STRING, 0, OPTIONAL, "Append Resampling internal info into this file") == XLAL_SUCCESS, XLAL_EFUNC );
XLALregBOOLUserStruct ( help, 'h', UVAR_HELP, "Print this message");
XLALregSTRINGUserStruct ( FstatMethod, 0, UVAR_OPTIONAL, XLALFstatMethodHelpString() );
XLALregREALUserStruct ( Freq, 0, UVAR_OPTIONAL, "Search frequency in Hz" );
XLALregREALUserStruct ( f1dot, 0, UVAR_OPTIONAL, "Search spindown f1dot in Hz/s" );
XLALregREALUserStruct ( FreqResolution, 0, UVAR_OPTIONAL, "Frequency resolution factor 'r' such that dFreq = 1/(r*T)" );
XLALregREALUserStruct ( Tseg, 0, UVAR_OPTIONAL, "Coherent segment length" );
XLALregINTUserStruct ( numSegments, 0, UVAR_OPTIONAL, "number of frequency bins to search" );
XLALregINTUserStruct ( numFreqBins, 0, UVAR_OPTIONAL, "number of frequency bins to search" );
XLALregREALUserStruct ( Tsft, 0, UVAR_DEVELOPER, "SFT length" );
XLAL_CHECK ( XLALRegisterUvarMember ( Tsft, REAL8, 0, DEVELOPER, "SFT length" ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALUserVarReadAllInput(argc, argv) == XLAL_SUCCESS, XLAL_EFUNC );
if (uvar->help) { // if help was requested, we're done here
......@@ -86,16 +95,12 @@ main ( int argc, char *argv[] )
XLAL_CHECK ( uvar->Tseg > uvar->Tsft, XLAL_EINVAL );
XLAL_CHECK ( uvar->Tsft > 1, XLAL_EINVAL );
XLAL_CHECK ( uvar->numFreqBins >= 1, XLAL_EINVAL );
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 );
fprintf ( stderr, "Tseg = %.1f d, numSegments = %" LAL_INT4_FORMAT ", Freq = %.1f Hz, f1dot = %.1e Hz/s, FreqResolution r = %f, numFreqBins = %" LAL_INT4_FORMAT " [dFreq = %.2e Hz, FreqBand = %.2e Hz]\n",
uvar->Tseg / 86400.0, uvar->numSegments, uvar->Freq, uvar->f1dot, uvar->FreqResolution, uvar->numFreqBins, dFreq, FreqBand );
// ---------- end: handle user input ----------
EphemerisData *ephem;
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 );
REAL8 memBase = XLALGetPeakHeapUsageMB();
REAL8 memEphem = memBase - maxMem0;
XLALPrintInfo ("mem(ephemeris) = %.1f MB\n", memEphem );
// ----- setup injection and data parameters
LALStringVector *detNames = NULL;
......@@ -144,6 +149,13 @@ main ( int argc, char *argv[] )
Doppler.ecc = ecc;
Doppler.asini = asini;
FILE *fpInfo = NULL;
if ( FstatMethod == FMETHOD_RESAMP_GENERIC )
{
XLAL_CHECK ( (fpInfo = fopen (uvar->outputInfo, "ab")) != NULL, XLAL_ESYS, "Failed to open '%s' for appending\n", uvar->outputInfo );
XLALAppendResampInfo2File ( fpInfo, NULL ); // create header comment line
}
// ----- setup optional Fstat arguments
FstatOptionalArgs optionalArgs = FstatOptionalArgsDefaults;
......@@ -160,12 +172,57 @@ main ( int argc, char *argv[] )
XLAL_CHECK ( (inputs = XLALCreateFstatInputVector ( uvar->numSegments )) != NULL, XLAL_EFUNC );
for ( INT4 l = 0; l < uvar->numSegments; l ++ )
{
XLAL_CHECK ( (inputs->data[l] = XLALCreateFstatInput ( catalogs[l], minCoverFreq, maxCoverFreq, dFreq, ephem, &optionalArgs )) != NULL, XLAL_EFUNC );
if ( l == 0 ) {
sharedWorkspace = XLALGetSharedFstatWorkspace ( inputs->data[0] );
}
optionalArgs.sharedWorkspace = sharedWorkspace;
}
// randomize numFreqBins
UINT4 numFreqBins_i = numFreqBinsMin + (UINT4)round ( 1.0 * (numFreqBinsMax - numFreqBinsMin) * rand() / RAND_MAX );
// randomize FreqResolution
REAL8 FreqResolution_i = FreqResolutionMin + 1.0 * ( FreqResolutionMax - FreqResolutionMin ) * rand() / RAND_MAX;
XLAL_CHECK ( (inputs = XLALCreateFstatInputVector ( uvar->numSegments )) != NULL, XLAL_EFUNC );
REAL8 dFreq = 1.0 / ( FreqResolution_i * uvar->Tseg );
REAL8 FreqBand = numFreqBins_i * dFreq;
fprintf ( stderr, "trial %d/%d: Tseg = %.1f d, numSegments = %d, Freq = %.1f Hz, f1dot = %.1e Hz/s, FreqResolution r = %f, numFreqBins = %d [dFreq = %.2e Hz, FreqBand = %.2e Hz]\n",
i+1, uvar->numTrials, uvar->Tseg / 86400.0, uvar->numSegments, uvar->Freq, uvar->f1dot, FreqResolution_i, numFreqBins_i, dFreq, FreqBand );
spinRange.fkdotBand[0] = FreqBand;
XLAL_CHECK ( XLALCWSignalCoveringBand ( &minCoverFreq, &maxCoverFreq, &startTime, &endTime, &spinRange, asini, Period, ecc ) == XLAL_SUCCESS, XLAL_EFUNC );
UINT4 numBinsSFT = ceil ( (maxCoverFreq - minCoverFreq) * uvar->Tsft + 2 * 8 );
REAL8 memSFTs = uvar->numSegments * numSFTsPerSeg * ( sizeof(SFTtype) + numBinsSFT * sizeof(COMPLEX8)) / 1e6;
// create per-segment input structs
for ( INT4 l = 0; l < uvar->numSegments; l ++ )
{
XLAL_CHECK ( (inputs->data[l] = XLALCreateFstatInput ( catalogs[l], minCoverFreq, maxCoverFreq, dFreq, ephem, &optionalArgs )) != NULL, XLAL_EFUNC );
if ( l == 0 ) {
sharedWorkspace = XLALGetSharedFstatWorkspace ( inputs->data[0] );
}
optionalArgs.sharedWorkspace = sharedWorkspace;
}
// ----- compute Fstatistics over segments
REAL8 tauFSumUnbuffered = 0;
REAL8 tauFSumBuffered = 0;
for ( INT4 l = 0; l < uvar->numSegments; l ++ )
{
REAL8 tic = XLALGetCPUTime();
XLAL_CHECK ( XLALComputeFstat ( &results, inputs->data[l], &Doppler, numFreqBins_i, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
REAL8 toc = XLALGetCPUTime();
tauFSumUnbuffered += ( toc - tic );
// ----- output more details if requested [only from first segment]
if ( (l == 0) && (fpInfo != NULL) ) {
XLALAppendResampInfo2File ( fpInfo, inputs->data[0] );
}
if ( uvar->runBuffered )
{
tic = XLALGetCPUTime();
XLAL_CHECK ( XLALComputeFstat ( &results, inputs->data[l], &Doppler, numFreqBins_i, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
toc = XLALGetCPUTime();
tauFSumBuffered += ( toc - tic );
}
// ----- compute Fstatistics over segments
FstatQuantities whatToCompute = (FSTATQ_2F | FSTATQ_2F_PER_DET);
......@@ -180,11 +237,20 @@ main ( int argc, char *argv[] )
XLAL_CHECK ( XLALComputeFstat ( &results[l], inputs->data[l], &Doppler, uvar->numFreqBins, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
REAL8 toc = XLALGetTimeOfDay();
tauFSumUnbuffered += ( toc - tic );
// ----- output more details if requested [first unbuffered segment]
if ( (l == 0) && (fpInfo != NULL) ) {
XLALAppendResampInfo2File ( fpInfo, inputs->data[0] );
}
// now call it with full buffering to get converged runtime per template (assuming many templates per skypoint or per binary params)
tic = XLALGetTimeOfDay();
XLAL_CHECK ( XLALComputeFstat ( &results[l], inputs->data[l], &Doppler, uvar->numFreqBins, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
toc = XLALGetTimeOfDay();
tauFSumBuffered += (toc - tic);
// ----- output more details if requested [first buffered segment]
if ( (l == 0) && (fpInfo != NULL) ) {
XLALAppendResampInfo2File ( fpInfo, inputs->data[0] );
}
} // for l < numSegments
REAL8 tauF1Buffered = tauFSumBuffered / ( uvar->numSegments * uvar->numFreqBins * numDetectors );
REAL8 tauF1Unbuffered = tauFSumUnbuffered / ( uvar->numSegments * uvar->numFreqBins * numDetectors );
......@@ -194,6 +260,9 @@ main ( int argc, char *argv[] )
XLALGetFstatMethodName ( FstatMethod ), tauF1Buffered, numDetectors * tauF1Buffered / numSFTsPerSeg, tauF1Unbuffered, memSFTs, memMaxCompute );
// ----- free memory ----------
if ( fpInfo != NULL ) {
fclose ( fpInfo );
}
XLALDestroyFstatInputVector ( inputs );
for ( INT4 l = 0; l < uvar->numSegments; l ++ )
{
......
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