Commit 65679b16 authored by Karl Wette's avatar Karl Wette Committed by Reinhard Prix
Browse files

ComputeFstat: better encapsulate shared workspace handling

- Move workspace into common Fstat data
- Move workspace ownership logic into XLALCreateFstatInput()
- Methods that use a workspace must supply a destructor function
- User now just needs to pass a previous FstatInput to optionalArgs;
  the FstatInput that allocated the workspace also destroys it
- Refs #2001
Original: 4168a79a1fcff721685f948f3b41a9a80fda801e
parent 16fa6a84
......@@ -155,7 +155,6 @@ typedef struct {
UINT4 nSFTs; /**< total number of SFTs */
LALStringVector *detectorIDs; /**< vector of detector IDs */
REAL4 NSegmentsInvX[PULSAR_MAX_DETECTORS]; /**< effective inverse number of segments per detector (needed for correct averaging in single-IFO F calculation) */
FstatWorkspace *sharedWorkspace; /**< resampling Fstat workspace shared across segments */
} UsefulStageVariables;
......@@ -1742,7 +1741,6 @@ int MAIN( int argc, char *argv[]) {
LALFree( fnameFstatVec1 );
}
XLALDestroyFstatWorkspace ( usefulParams.sharedWorkspace );
XLALDestroyFstatInputVector(Fstat_in_vec);
XLALDestroyFstatResults(Fstat_res);
......@@ -2060,9 +2058,9 @@ void SetUpSFTs( LALStatus *status, /**< pointer to LALStatus structure */
XLALPrintError("%s: XLALCreateFstatInput() failed with errno=%d", __func__, xlalErrno);
ABORT ( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
}
// re-use shared Resampling workspace from first segment for all segments
// re-use shared workspace from first segment for all segments
if ( k == 0 ) {
in->sharedWorkspace = optionalArgs.sharedWorkspace = XLALGetSharedFstatWorkspace ( (*p_Fstat_in_vec)->data[0] );
optionalArgs.prevInput = (*p_Fstat_in_vec)->data[0];
}
if ( k == 0 ) {
LogPrintf (LOG_NORMAL, "FstatMethod used: '%s'\n", XLALGetFstatInputMethodName( (*p_Fstat_in_vec)->data[k] ) );
......
......@@ -49,8 +49,14 @@ typedef struct {
SSBprecision SSBprec; // Barycentric transformation precision
FstatMethodType FstatMethod; // Method to use for computing the F-statistic
REAL8 dFreq; // Requested spacing of \f$\mathcal{F}\f$-statistic frequency bins.
void *workspace; // F-statistic method workspace
} FstatInput_Common;
// Pointers to function pointers which perform method-specific operations
typedef struct {
void (*workspace_dtor)(void *); // Workspace destructor function
} FstatInput_MethodFuncs;
// Input data specific to F-statistic methods
typedef struct tagFstatInput_Demod FstatInput_Demod;
typedef struct tagFstatInput_Resamp FstatInput_Resamp;
......@@ -59,6 +65,8 @@ typedef struct tagFstatInput_Resamp FstatInput_Resamp;
struct tagFstatInput {
FstatMethodType method; // Method to use for computing the F-statistic
FstatInput_Common common; // Common input data
BOOLEAN own_workspace; // True if we own 'common.workspace', and therefore must destroy it
FstatInput_MethodFuncs method_funcs; // Function pointers for F-statistic method
FstatInput_Demod* demod; // Demodulation input data
FstatInput_Resamp* resamp; // Resampling input data
};
......@@ -110,7 +118,7 @@ const FstatOptionalArgs FstatOptionalArgsDefaults = {
.injectSources = NULL,
.injectSqrtSX = NULL,
.assumeSqrtSX = NULL,
.sharedWorkspace = NULL
.prevInput = NULL
};
// ---------- Include F-statistic method implementations ---------- //
......@@ -320,6 +328,16 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Cata
common->FstatMethod = optArgs->FstatMethod;
common->dFreq = dFreq;
// Determine whether we can re-used workspace from a previous call to XLALCreateFstatInput()
if ( optArgs->prevInput != NULL ) {
XLAL_CHECK_NULL( optArgs->prevInput->method == input->method, XLAL_EFAILED );
common->workspace = optArgs->prevInput->common.workspace;
input->own_workspace = 0; // We do not own this workspace, the original owner will destroy it
} else {
common->workspace = NULL;
input->own_workspace = 1; // We own this workspace, and therefore have the responsibility to destroy it
}
// Determine the time baseline of an SFT
const REAL8 Tsft = 1.0 / SFTcatalog->data[0].header.deltaF;
......@@ -451,19 +469,19 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Cata
// Call the appropriate method function to setup their input data structures
// - The method input data structures are expected to take ownership of the
// SFTs, which is why 'input->common' does not retain a pointer to them
FstatInput_MethodFuncs *funcs = &input->method_funcs;
if ( input->demod != NULL )
{
XLAL_CHECK_NULL ( SetupFstatInput_Demod ( input->demod, common, multiSFTs ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK_NULL ( SetupFstatInput_Demod ( input->demod, common, funcs, multiSFTs ) == XLAL_SUCCESS, XLAL_EFUNC );
}
else if ( input->resamp != NULL )
{
XLAL_CHECK_NULL ( SetupFstatInput_Resamp ( input->resamp, common, multiSFTs, optArgs->sharedWorkspace ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK_NULL ( SetupFstatInput_Resamp ( input->resamp, common, funcs, multiSFTs ) == XLAL_SUCCESS, XLAL_EFUNC );
}
else
{
XLAL_ERROR_NULL ( XLAL_EFAILED, "Invalid FstatInput struct passed to %s()", __func__ );
}
multiSFTs = NULL;
// Cleanup
XLALDestroyMultiPSDVector ( runningMedian );
......@@ -686,6 +704,11 @@ XLALDestroyFstatInput ( FstatInput* input ///< [in] \c FstatInput structur
XLALDestroyMultiNoiseWeights ( input->common.multiNoiseWeights );
XLALDestroyMultiDetectorStateSeries ( input->common.multiDetectorStates );
if ( input->own_workspace && input->common.workspace != NULL ) {
// Free method-specific workspace using destructor function
(input->method_funcs.workspace_dtor) ( input->common.workspace );
}
if (input->demod != NULL)
{
DestroyFstatInput_Demod ( input->demod );
......
......@@ -66,11 +66,6 @@ extern "C" {
///
typedef struct tagFstatInput FstatInput;
///
/// pre-allocated 'workspace' memory for Resampling Fstat computation; can be shared over segments via XLALTakeOverFstatWorkspace()
///
typedef struct tagFstatWorkspace FstatWorkspace;
///
/// A vector of XLALComputeFstat() input data structures, for e.g. computing the
/// \f$\mathcal{F}\f$-statistic for multiple segments.
......@@ -161,7 +156,7 @@ typedef struct tagFstatOptionalArgs {
PulsarParamsVector *injectSources; ///< Vector of parameters of CW signals to simulate and inject.
MultiNoiseFloor *injectSqrtSX; ///< Single-sided PSD values for fake Gaussian noise to be added to SFT data.
MultiNoiseFloor *assumeSqrtSX; ///< Single-sided PSD values to be used for computing SFT noise weights instead of from a running median of the SFTs themselves.
FstatWorkspace *sharedWorkspace; ///< use this shared workspace instead of creating our own one
FstatInput *prevInput; ///< An \c FstatInput structure from a previous call to XLALCreateFstatInput(); may contain common workspace data than can be re-used to save memory.
} FstatOptionalArgs;
#ifdef SWIG // SWIG interface directives
SWIGLAL(COPY_CONSTRUCTOR(tagFstatOptionalArgs));
......@@ -321,9 +316,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 );
void XLALDestroyFstatInput ( FstatInput* input );
void XLALDestroyFstatResults ( FstatResults* Fstats );
int XLALAdd4ToFstatResults ( FstatResults* Fstats );
......
......@@ -301,13 +301,15 @@ DestroyFstatInput_Demod ( FstatInput_Demod* demod )
static int
SetupFstatInput_Demod ( FstatInput_Demod *demod,
const FstatInput_Common *common,
FstatInput_Common *common,
FstatInput_MethodFuncs* funcs,
MultiSFTVector *multiSFTs
)
{
// Check input
XLAL_CHECK ( common != NULL, XLAL_EFAULT );
XLAL_CHECK ( demod != NULL, XLAL_EFAULT );
XLAL_CHECK ( common != NULL, XLAL_EFAULT );
XLAL_CHECK ( funcs != NULL, XLAL_EFAULT );
XLAL_CHECK ( multiSFTs != NULL, XLAL_EFAULT );
// Save pointer to SFTs
......
......@@ -72,7 +72,7 @@ typedef struct tagResampTimingInfo
REAL8 tauF1NoBuf; // Resampling timing 'constant': Fstat time per template per detector for an 'unbuffered' usage (different skypos and numFreqBins)
} ResampTimingInfo;
struct tagFstatWorkspace
typedef struct tagResampWorkspace
{
// intermediate quantities to interpolate and operate on SRC-frame timeseries
COMPLEX8Vector *TStmp1_SRC; // can hold a single-detector SRC-frame spindown-corrected timeseries [without zero-padding]
......@@ -94,7 +94,7 @@ struct tagFstatWorkspace
UINT4 numFreqBinsAlloc; // internal: keep track of allocated length of frequency-arrays
ResampTimingInfo timingInfo; // temporary storage for collecting timing data
};
} ResampWorkspace;
struct tagFstatInput_Resamp
{
......@@ -107,10 +107,6 @@ struct tagFstatInput_Resamp
MultiCOMPLEX8TimeSeries *multiTimeSeries_SRC_a; // multi-detector SRC-frame timeseries, multiplied by AM function a(t)
MultiCOMPLEX8TimeSeries *multiTimeSeries_SRC_b; // multi-detector SRC-frame timeseries, multiplied by AM function b(t)
// ----- workspace -----
FstatWorkspace *ws; // 'workspace': pre-allocated vectors used to store intermediate results
BOOLEAN ownThisWorkspace; // flag whether we 'own' or share this workspace (ie who is responsible for freeing it)
};
......@@ -129,32 +125,29 @@ XLALBarycentricResampleMultiCOMPLEX8TimeSeries ( FstatInput_Resamp *resamp,
);
static int
XLALComputeFaFb_Resamp ( FstatWorkspace *ws,
XLALComputeFaFb_Resamp ( ResampWorkspace *ws,
const PulsarDopplerParams thisPoint,
REAL8 dFreq,
const COMPLEX8TimeSeries *TimeSeries_SRC_a,
const COMPLEX8TimeSeries *TimeSeries_SRC_b
);
static FstatWorkspace *XLALCreateFstatWorkspace ( UINT4 numSamplesSRC, UINT4 numSamplesFFT );
// pseudo-internal: don't export API but allow using them from test/benchmark codes
int XLALAppendResampInfo2File ( FILE *fp, const FstatInput *input );
int XLALGetResampTimingInfo ( REAL8 *tauF1NoBuf, REAL8 *tauF1Buf, const FstatInput *input );
// ==================== function definitions ====================
// ---------- exported API functions ----------
///
/// Create a new workspace with given time samples in SRC frame 'numSamplesSRC' (holds time-series for spindown-correction)
/// and given total number of time-samples for FFTing (includes zero-padding for frequency-resolution)
///
static FstatWorkspace *
XLALCreateFstatWorkspace ( UINT4 numSamplesSRC,
static ResampWorkspace *
XLALCreateResampWorkspace ( UINT4 numSamplesSRC,
UINT4 numSamplesFFT
)
{
FstatWorkspace *ws;
ResampWorkspace *ws;
XLAL_CHECK_NULL ( (ws = XLALCalloc ( 1, sizeof(*ws))) != NULL, XLAL_ENOMEM );
XLAL_CHECK_NULL ( (ws->TStmp1_SRC = XLALCreateCOMPLEX8Vector ( numSamplesSRC )) != NULL, XLAL_EFUNC );
......@@ -170,35 +163,12 @@ XLALCreateFstatWorkspace ( UINT4 numSamplesSRC,
ws->numSamplesFFT = numSamplesFFT;
return ws;
} // XLALCreateFstatWorkspace()
///
/// Function to extract a workspace from a resampling setup, which can be passed in FstatOptionalArgs to be shared by various setups
/// in order to save memory. Note, when using this, you need to free this workspace yourself at the end using XLALDestroyFstatWorkspace().
/// Note: Demod methods don't use a workspace, so NULL (without error) is returned in this case.
///
FstatWorkspace *
XLALGetSharedFstatWorkspace ( FstatInput *input //!< [in,out] Fstat input structure to extract shared workspace from
)
{
XLAL_CHECK_NULL ( input != NULL, XLAL_EINVAL );
} // XLALCreateResampWorkspace()
if ( input->resamp == NULL ) {
return NULL;
}
input->resamp->ownThisWorkspace = 0; // the caller now owns the workspace and has to free it
return input->resamp->ws;
} // XLALGetSharedFstatWorkspace()
void
XLALDestroyFstatWorkspace ( FstatWorkspace *ws )
static void
XLALDestroyResampWorkspace ( void *workspace )
{
if ( ws == NULL ) {
return;
}
ResampWorkspace *ws = (ResampWorkspace*) workspace;
XLALDestroyCOMPLEX8Vector ( ws->TStmp1_SRC );
XLALDestroyCOMPLEX8Vector ( ws->TStmp2_SRC );
......@@ -219,16 +189,17 @@ XLALDestroyFstatWorkspace ( FstatWorkspace *ws )
XLALFree ( ws );
return;
} // XLALDestroyFstatWorkspace()
} // XLALDestroyResampWorkspace()
// return resampling timing coefficients 'tauF1' for buffered and unbuffered calls
int
XLALGetResampTimingInfo ( REAL8 *tauF1NoBuf, REAL8 *tauF1Buf, const FstatInput *input )
{
XLAL_CHECK ( input != NULL, XLAL_EINVAL );
XLAL_CHECK ( input->resamp != NULL, XLAL_EINVAL );
XLAL_CHECK ( input->resamp->ws != NULL, XLAL_EINVAL );
const ResampTimingInfo *ti = &(input->resamp->ws->timingInfo);
XLAL_CHECK ( input->method >= FMETHOD_RESAMP_GENERIC, XLAL_EINVAL );
const FstatInput_Common *common = &input->common;
const ResampWorkspace *ws = (ResampWorkspace*) common->workspace;
const ResampTimingInfo *ti = &(ws->timingInfo);
(*tauF1NoBuf) = ti->tauF1NoBuf;
(*tauF1Buf) = ti->tauF1Buf;
......@@ -251,10 +222,10 @@ XLALAppendResampInfo2File ( FILE *fp, const FstatInput *input )
"tauTotal", "tauFFT", "tauBary", "tauSpin", "tauAM", "tauNorm", "tauFab2F", "tauMem", "tauSumFabX", "tauF1NoBuf", "tauF1Buf" );
return XLAL_SUCCESS;
}
XLAL_CHECK ( input->resamp != NULL, XLAL_EINVAL );
XLAL_CHECK ( input->method >= FMETHOD_RESAMP_GENERIC, XLAL_EINVAL );
const FstatInput_Resamp *resamp = input->resamp;
const FstatInput_Common *common = &input->common;
const FstatWorkspace *ws = resamp->ws;
const ResampWorkspace *ws = (ResampWorkspace*) common->workspace;
fprintf (fp, "%10d %10d", ws->numFreqBinsOut, ws->numSamplesFFT );
UINT4 numSamples_DETX0 = resamp->multiTimeSeries_DET->data[0]->data->length;
......@@ -281,11 +252,6 @@ DestroyFstatInput_Resamp ( FstatInput_Resamp* resamp )
XLALDestroyMultiCOMPLEX8TimeSeries ( resamp->multiTimeSeries_SRC_a );
XLALDestroyMultiCOMPLEX8TimeSeries ( resamp->multiTimeSeries_SRC_b );
// ----- free workspace
if ( resamp->ownThisWorkspace ) {
XLALDestroyFstatWorkspace ( resamp->ws );
}
XLALFree ( resamp );
return;
......@@ -294,15 +260,16 @@ DestroyFstatInput_Resamp ( FstatInput_Resamp* resamp )
static int
SetupFstatInput_Resamp ( FstatInput_Resamp *resamp,
const FstatInput_Common *common,
MultiSFTVector *multiSFTs,
FstatWorkspace *sharedWorkspace
)
FstatInput_Common *common,
FstatInput_MethodFuncs* funcs,
MultiSFTVector *multiSFTs
)
{
// Check input
XLAL_CHECK(common != NULL, XLAL_EFAULT);
XLAL_CHECK(resamp != NULL, XLAL_EFAULT);
XLAL_CHECK(multiSFTs != NULL, XLAL_EFAULT);
XLAL_CHECK ( resamp != NULL, XLAL_EFAULT );
XLAL_CHECK ( common != NULL, XLAL_EFAULT );
XLAL_CHECK ( funcs != NULL, XLAL_EFAULT );
XLAL_CHECK ( multiSFTs != NULL, XLAL_EFAULT );
// Convert SFTs into heterodyned complex timeseries [in detector frame]
XLAL_CHECK ( (resamp->multiTimeSeries_DET = XLALMultiSFTVectorToCOMPLEX8TimeSeries ( multiSFTs )) != NULL, XLAL_EFUNC );
......@@ -358,8 +325,9 @@ SetupFstatInput_Resamp ( FstatInput_Resamp *resamp,
} // for X < numDetectors
// ---- re-use shared workspace, or allocate here ----------
if ( sharedWorkspace != NULL )
if ( common->workspace != NULL )
{
ResampWorkspace *sharedWorkspace = (ResampWorkspace*) common->workspace;
XLAL_CHECK ( numSamplesFFT == sharedWorkspace->numSamplesFFT, XLAL_EINVAL, "Shared workspace of different frequency resolution: numSamplesFFT = %d != %d\n",
sharedWorkspace->numSamplesFFT, numSamplesFFT );
......@@ -372,13 +340,11 @@ SetupFstatInput_Resamp ( FstatInput_Resamp *resamp,
XLAL_CHECK ( (sharedWorkspace->SRCtimes_DET->data = XLALRealloc ( sharedWorkspace->SRCtimes_DET->data, numSamplesMax_SRC * sizeof(REAL8) )) != NULL, XLAL_ENOMEM );
sharedWorkspace->SRCtimes_DET->length = numSamplesMax_SRC;
}
resamp->ws = sharedWorkspace;
resamp->ownThisWorkspace = 0;
} // end: if shared workspace given
else
{
XLAL_CHECK ( ( resamp->ws = XLALCreateFstatWorkspace ( numSamplesMax_SRC, numSamplesFFT )) != NULL, XLAL_EFUNC );
resamp->ownThisWorkspace = 1;
XLAL_CHECK ( ( common->workspace = XLALCreateResampWorkspace ( numSamplesMax_SRC, numSamplesFFT )) != NULL, XLAL_EFUNC );
funcs->workspace_dtor = XLALDestroyResampWorkspace;
} // end: if we create our own workspace
return XLAL_SUCCESS;
......@@ -408,10 +374,12 @@ 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" );
ResampWorkspace *ws = (ResampWorkspace*) common->workspace;
#ifdef COLLECT_TIMING
// collect internal timing info
XLAL_INIT_MEM ( resamp->ws->timingInfo );
ResampTimingInfo *ti = &(resamp->ws->timingInfo);
XLAL_INIT_MEM ( ws->timingInfo );
ResampTimingInfo *ti = &(ws->timingInfo);
REAL8 ticStart,tocEnd;
ticStart = XLALGetCPUTime();
REAL8 tic,toc;
......@@ -432,8 +400,6 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
(GPSDIFF( resamp->prev_doppler.tp, thisPoint.tp ) == 0 ) &&
(resamp->prev_doppler.argp == thisPoint.argp);
FstatWorkspace *ws = resamp->ws;
// ----- not same skypos+binary+refTime? --> re-compute SRC-frame timeseries, AM-coeffs and store in buffer
#ifdef COLLECT_TIMING
tic = XLALGetCPUTime();
......@@ -515,8 +481,8 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
const COMPLEX8TimeSeries *TimeSeriesX_SRC_a = multiTimeSeries_SRC_a->data[X];
const COMPLEX8TimeSeries *TimeSeriesX_SRC_b = multiTimeSeries_SRC_b->data[X];
// 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_a, TimeSeriesX_SRC_b ) == XLAL_SUCCESS, XLAL_EFUNC );
// compute {Fa^X(f_k), Fb^X(f_k)}: results returned via workspace ws
XLAL_CHECK ( XLALComputeFaFb_Resamp ( ws, thisPoint, common->dFreq, TimeSeriesX_SRC_a, TimeSeriesX_SRC_b ) == XLAL_SUCCESS, XLAL_EFUNC );
#ifdef COLLECT_TIMING
tic = XLALGetCPUTime();
......@@ -628,7 +594,7 @@ ComputeFstat_Resamp ( FstatResults* Fstats,
static int
XLALComputeFaFb_Resamp ( FstatWorkspace *restrict ws, //!< [in,out] pre-allocated 'workspace' for temporary and output quantities
XLALComputeFaFb_Resamp ( ResampWorkspace *restrict ws, //!< [in,out] pre-allocated 'workspace' for temporary and output quantities
const PulsarDopplerParams thisPoint, //!< [in] Doppler point to compute {FaX,FbX} for
REAL8 dFreq, //!< [in] output frequency resolution
const COMPLEX8TimeSeries * restrict TimeSeries_SRC_a, //!< [in] SRC-frame single-IFO timeseries * a(t)
......@@ -797,6 +763,8 @@ XLALBarycentricResampleMultiCOMPLEX8TimeSeries ( FstatInput_Resamp *resamp, //
XLAL_CHECK ( resamp->multiTimeSeries_SRC_a != NULL, XLAL_EINVAL );
XLAL_CHECK ( resamp->multiTimeSeries_SRC_b != NULL, XLAL_EINVAL );
ResampWorkspace *ws = (ResampWorkspace*) common->workspace;
UINT4 numDetectors = resamp->multiTimeSeries_DET->length;
XLAL_CHECK ( resamp->multiTimeSeries_SRC_a->length == numDetectors, XLAL_EINVAL, "Inconsistent number of detectors tsDET(%d) != tsSRC(%d)\n", numDetectors, resamp->multiTimeSeries_SRC_a->length );
XLAL_CHECK ( resamp->multiTimeSeries_SRC_b->length == numDetectors, XLAL_EINVAL, "Inconsistent number of detectors tsDET(%d) != tsSRC(%d)\n", numDetectors, resamp->multiTimeSeries_SRC_b->length );
......@@ -843,7 +811,7 @@ XLALBarycentricResampleMultiCOMPLEX8TimeSeries ( FstatInput_Resamp *resamp, //
// shorthand pointers: output
COMPLEX8TimeSeries *TimeSeries_SRCX_a = resamp->multiTimeSeries_SRC_a->data[X];
COMPLEX8TimeSeries *TimeSeries_SRCX_b = resamp->multiTimeSeries_SRC_b->data[X];
REAL8Vector *ti_DET = resamp->ws->SRCtimes_DET;
REAL8Vector *ti_DET = ws->SRCtimes_DET;
// useful shorthands
REAL8 refTime8 = GPSGETREAL8 ( &SRCtimesX->refTime );
......@@ -873,7 +841,7 @@ XLALBarycentricResampleMultiCOMPLEX8TimeSeries ( FstatInput_Resamp *resamp, //
memset ( TimeSeries_SRCX_a->data->data, 0, TimeSeries_SRCX_a->data->length * sizeof(TimeSeries_SRCX_a->data->data[0]) );
memset ( TimeSeries_SRCX_b->data->data, 0, TimeSeries_SRCX_b->data->length * sizeof(TimeSeries_SRCX_b->data->data[0]) );
// make sure detector-frame timesteps to interpolate to are initialized to 0, in case of gaps
memset ( resamp->ws->SRCtimes_DET->data, 0, resamp->ws->SRCtimes_DET->length * sizeof(resamp->ws->SRCtimes_DET->data[0]) );
memset ( ws->SRCtimes_DET->data, 0, ws->SRCtimes_DET->length * sizeof(ws->SRCtimes_DET->data[0]) );
REAL8 tStart_DET_0 = GPSGETREAL8 ( &(Timestamps_DETX->data[0]) );// START time of the SFT at the detector
......@@ -922,8 +890,8 @@ XLALBarycentricResampleMultiCOMPLEX8TimeSeries ( FstatInput_Resamp *resamp, //
// apply AM coefficients a(t), b(t) to SRC frame timeseries [alternate sign to get final FFT return DC in the middle]
REAL4 signum = signumLUT [ (iSRC_al_j % 2) ]; // alternating sign, avoid branching
ei2piphase *= signum;
resamp->ws->TStmp1_SRC->data [ iSRC_al_j ] = ei2piphase * a_al;
resamp->ws->TStmp2_SRC->data [ iSRC_al_j ] = ei2piphase * b_al;
ws->TStmp1_SRC->data [ iSRC_al_j ] = ei2piphase * a_al;
ws->TStmp2_SRC->data [ iSRC_al_j ] = ei2piphase * b_al;
} // for j < numSamples_SRC_al
} // for alpha < numSFTsX
......@@ -938,8 +906,8 @@ XLALBarycentricResampleMultiCOMPLEX8TimeSeries ( FstatInput_Resamp *resamp, //
// apply heterodyne correction and AM-functions a(t) and b(t) to interpolated timeseries
for ( UINT4 j = 0; j < numSamples_SRCX; j ++ )
{
TimeSeries_SRCX_b->data->data[j] = TimeSeries_SRCX_a->data->data[j] * resamp->ws->TStmp2_SRC->data[j];
TimeSeries_SRCX_a->data->data[j] *= resamp->ws->TStmp1_SRC->data[j];
TimeSeries_SRCX_b->data->data[j] = TimeSeries_SRCX_a->data->data[j] * ws->TStmp2_SRC->data[j];
TimeSeries_SRCX_a->data->data[j] *= ws->TStmp1_SRC->data[j];
} // for j < numSamples_SRCX
} // for X < numDetectors
......
......@@ -195,7 +195,6 @@ main ( int argc, char *argv[] )
optionalArgs.injectSqrtSX = &injectSqrtSX;
optionalArgs.FstatMethod = FstatMethod;
FstatWorkspace *sharedWorkspace = NULL;
FstatInputVector *inputs;
FstatQuantities whatToCompute = (FSTATQ_2F | FSTATQ_2F_PER_DET);
FstatResults *results = NULL;
......@@ -228,9 +227,8 @@ main ( int argc, char *argv[] )
{
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.prevInput = inputs->data[0];
}
optionalArgs.sharedWorkspace = sharedWorkspace;
}
// ----- compute Fstatistics over segments
......@@ -261,8 +259,6 @@ main ( int argc, char *argv[] )
XLALGetFstatInputMethodName ( inputs->data[0] ), tauF1Buf_i, tauF1NoBuf_i, memSFTs, memMaxCompute );
XLALDestroyFstatInputVector ( inputs );
XLALDestroyFstatWorkspace ( sharedWorkspace );
optionalArgs.sharedWorkspace = NULL;
} // for i < numTrials
......
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