31 #include <lal/FFTWMutex.h>
32 #include <lal/Factorial.h>
33 #include <lal/LFTandTSutils.h>
34 #include <lal/LogPrintf.h>
35 #include <lal/SinCosLUT.h>
36 #include <lal/TimeSeries.h>
37 #include <lal/Units.h>
57 typedef struct tagResampGenericWorkspace {
156 fftwf_destroy_plan( resamp->
fftplan );
217 REAL8 TspanX = numSamples_DETX * dt_DET;
218 TspanXMax =
fmax( TspanXMax, TspanX );
220 UINT4 decimateFFT = (
UINT4 )ceil( TspanXMax / TspanFFT );
221 if ( decimateFFT > 1 ) {
222 XLALPrintWarning(
"WARNING: Frequency spacing larger than 1/Tspan, we'll internally decimate FFT frequency bins by a factor of %" LAL_UINT4_FORMAT "\n", decimateFFT );
224 TspanFFT *= decimateFFT;
227 UINT4 numSamplesFFT0 = (
UINT4 ) ceil( TspanFFT / dt_DET );
228 UINT4 numSamplesFFT = 0;
230 numSamplesFFT = (
UINT4 ) pow( 2, ceil( log2( numSamplesFFT0 ) ) );
232 numSamplesFFT = (
UINT4 ) 2 * ceil( numSamplesFFT0 / 2 );
235 REAL8 dt_SRC = TspanFFT / numSamplesFFT;
250 UINT4 numSamplesMax_SRC = 0;
254 XLAL_CHECK( dt_DET == dt_DETX,
XLAL_EINVAL,
"Input timeseries must have identical 'deltaT(X=%d)' (%.16g != %.16g)\n", X, dt_DET, dt_DETX );
257 XLAL_CHECK( fabs( fHet - fHetX ) <
LAL_REAL8_EPS * fHet,
XLAL_EINVAL,
"Input timeseries must have identical heterodyning frequency 'f0(X=%d)' (%.16g != %.16g)\n", X, fHet, fHetX );
265 UINT4 numSamples_SRCX = (
UINT4 )ceil( numSamples_DETX * dt_DET / dt_SRC );
270 numSamplesMax_SRC =
MYMAX( numSamplesMax_SRC, numSamples_SRCX );
273 XLAL_CHECK( numSamplesFFT >= numSamplesMax_SRC,
XLAL_EFAILED,
"[numSamplesFFT = %d] < [numSamplesMax_SRC = %d]\n", numSamplesFFT, numSamplesMax_SRC );
312 int fft_plan_flags = FFTW_MEASURE;
313 double fft_plan_timeout = FFTW_NO_TIMELIMIT ;
314 char *wisdom_filename;
315 static int tried_wisdom = 0;
319 wisdom_filename = getenv(
"FFTWF_WISDOM_FILENAME" );
320 if ( wisdom_filename && !tried_wisdom ) {
321 FILE *
fp = fopen( wisdom_filename,
"r" );
323 XLALPrintWarning(
"WARNING: Couldn't open wisdom file '%s'\n", wisdom_filename );
324 }
else if ( fftwf_import_wisdom_from_file(
fp ) ) {
325 XLALPrintInfo(
"INFO: imported wisdom from file '%s'\n", wisdom_filename );
328 XLALPrintWarning(
"WARNING: Couldn't import wisdom from file '%s'\n", wisdom_filename );
334 fftw_set_timelimit( fft_plan_timeout );
385 REAL8 ticStart = 0, tocEnd = 0;
386 REAL8 tic = 0, toc = 0;
387 if ( collectTiming ) {
406 if ( collectTiming ) {
442 if ( collectTiming ) {
444 Tau->Mem = ( toc - tic );
461 if ( collectTiming ) {
466 for (
UINT4 k = 0;
k < numFreqBins;
k++ ) {
473 for (
UINT4 k = 0;
k < numFreqBins;
k++ ) {
479 if ( collectTiming ) {
481 Tau->SumFabX += ( toc - tic );
492 for (
UINT4 k = 0;
k < numFreqBins;
k ++ ) {
497 if ( collectTiming ) {
499 Tau->Fab2F += ( toc - tic );
504 if ( collectTiming ) {
516 for (
UINT4 k = 0;
k < numFreqBins;
k++ ) {
524 if ( collectTiming ) {
526 Tau->Fab2F += ( toc - tic );
555 if ( collectTiming ) {
562 Tau->Total = ( tocEnd - ticStart );
570 REAL8 Tau_buffer = Tau->Bary;
573 REAL8 tauF_eff = Tau->Total / NFbin;
574 REAL8 tauF_core = ( Tau->Total - Tau_buffer ) / NFbin;
577 REAL8 tau0_Fbin = ( Tau->Copy + Tau->Norm + Tau->SumFabX + Tau->Fab2F ) / NFbin;
583 #define updateAvgF(q) tiGen->q = ((tiGen->q *(tiGen->NCalls-1) + q)/(tiGen->NCalls))
589 #define updateAvgRS(q) tiRS->q = ((tiRS->q *(tiGen->NCalls-1) + q)/(tiGen->NCalls))
595 if ( Tau->BufferRecomputed ) {
597 REAL8 tauF_buffer = Tau_buffer / NFbin;
620 XLAL_CHECK( ( resamp != NULL ) && ( ws != NULL ) && ( TimeSeries_SRC_a != NULL ) && ( TimeSeries_SRC_b != NULL ),
XLAL_EINVAL );
627 REAL8 fHet = TimeSeries_SRC_a->f0;
628 REAL8 dt_SRC = TimeSeries_SRC_a->deltaT;
631 REAL8 freqShift = remainder( FreqOut0 - fHet, dFreq );
633 XLAL_CHECK( FreqOut0 >= fMinFFT,
XLAL_EDOM,
"Lowest output frequency outside the available frequency band: [FreqOut0 = %.16g] < [fMinFFT = %.16g]\n", FreqOut0, fMinFFT );
634 UINT4 offset_bins = (
UINT4 ) lround( ( FreqOut0 - fMinFFT ) / dFreqFFT );
635 UINT4 maxOutputBin = offset_bins + ( numFreqBins - 1 ) * resamp->
decimateFFT;
636 XLAL_CHECK( maxOutputBin < resamp->numSamplesFFT,
XLAL_EDOM,
"Highest output frequency bin outside available band: [maxOutputBin = %d] >= [numSamplesFFT = %d]\n", maxOutputBin, resamp->
numSamplesFFT );
640 REAL8 tic = 0, toc = 0;
645 if ( collectTiming ) {
653 if ( collectTiming ) {
655 tiRS->
Tau.
Spin += ( toc - tic );
662 if ( collectTiming ) {
664 tiRS->
Tau.
FFT += ( toc - tic );
668 for (
UINT4 k = 0;
k < numFreqBins;
k++ ) {
672 if ( collectTiming ) {
674 tiRS->
Tau.
Copy += ( toc - tic );
682 if ( collectTiming ) {
684 tiRS->
Tau.
Spin += ( toc - tic );
691 if ( collectTiming ) {
693 tiRS->
Tau.
FFT += ( toc - tic );
697 for (
UINT4 k = 0;
k < numFreqBins;
k++ ) {
701 if ( collectTiming ) {
703 tiRS->
Tau.
Copy += ( toc - tic );
709 for (
UINT4 k = 0;
k < numFreqBins;
k++ ) {
710 REAL8 f_k = FreqOut0 +
k * dFreq;
711 REAL8 cycles = - f_k * dtauX;
712 REAL4 sinphase, cosphase;
719 if ( collectTiming ) {
721 tiRS->
Tau.
Norm += ( toc - tic );
743 while ( ( s_max > 0 ) && ( doppler->fkdot[s_max] == 0 ) ) {
748 UINT4 numSamplesIn = xIn->data->length;
754 for (
UINT4 j = 0;
j < numSamplesIn;
j ++ ) {
756 REAL8 Dtau_alpha_j = Dtau0 + taup_j;
758 REAL8 cycles = - freqShift * taup_j;
760 REAL8 Dtau_pow_kp1 = Dtau_alpha_j;
761 for (
UINT4 k = 1;
k <= s_max;
k++ ) {
762 Dtau_pow_kp1 *= Dtau_alpha_j;
763 cycles += -
LAL_FACT_INV[
k + 1] * doppler->fkdot[
k] * Dtau_pow_kp1;
766 REAL4 cosphase, sinphase;
771 xOut[
j] = em2piphase * xIn->data->data[
j];
818 REAL8 tic = 0, toc = 0;
822 if ( same_skypos && same_refTime && same_binary ) {
829 if ( collectTiming ) {
836 if ( !( same_skypos && same_refTime ) ) {
858 if ( thisPoint->
asini > 0 ) {
873 const REAL4 signumLUT[2] = {1, -1};
898 XLAL_CHECK( numSamples_DETX > 0,
XLAL_EINVAL,
"Input timeseries for detector X=%d has zero samples. Can't handle that!\n", X );
901 XLAL_CHECK( fabs( fHet - fHetX ) <
LAL_REAL8_EPS * fHet,
XLAL_EINVAL,
"Input timeseries must have identical heterodyning frequency 'f0(X=%d)' (%.16g != %.16g)\n", X, fHet, fHetX );
905 TimeSeries_SRCX_a->
f0 = fHet;
906 TimeSeries_SRCX_b->
f0 = fHet;
930 REAL8 tStart_SRC_al = tMid_SRC_al - 0.5 *
Tsft * Tdot_al;
931 REAL8 tEnd_SRC_al = tMid_SRC_al + 0.5 *
Tsft * Tdot_al;
934 REAL8 tMid_DET_al = tStart_DET_al + 0.5 *
Tsft;
937 UINT4 iStart_SRC_al = lround( ( tStart_SRC_al - tStart_SRC_0 ) / dt_SRC );
938 UINT4 iEnd_SRC_al = lround( ( tEnd_SRC_al - tStart_SRC_0 ) / dt_SRC );
941 iStart_SRC_al =
MYMIN( iStart_SRC_al, numSamples_SRCX - 1 );
942 iEnd_SRC_al =
MYMIN( iEnd_SRC_al, numSamples_SRCX - 1 );
943 UINT4 numSamplesSFT_SRC_al = iEnd_SRC_al - iStart_SRC_al + 1;
947 for (
UINT4 j = 0;
j < numSamplesSFT_SRC_al;
j++ ) {
948 UINT4 iSRC_al_j = iStart_SRC_al +
j;
952 REAL8 t_SRC = tStart_SRC_0 + iSRC_al_j * dt_SRC;
953 ti_DET->
data [ iSRC_al_j ] = tMid_DET_al + ( t_SRC - tMid_SRC_al ) / Tdot_al;
956 REAL8 tDiff = iSRC_al_j * dt_SRC + ( tStart_DET_0 - ti_DET->
data [ iSRC_al_j ] );
957 REAL8 cycles = fmod( fHet * tDiff, 1.0 );
960 REAL4 cosphase, sinphase;
965 REAL4 signum = signumLUT [( iSRC_al_j % 2 ) ];
966 ei2piphase *= signum;
977 ti_DET->
length = bak_length;
980 for (
UINT4 j = 0;
j < numSamples_SRCX;
j ++ ) {
987 if ( collectTiming ) {
989 Tau->
Bary = ( toc - tic );
998 double *planGenTimeoutSeconds
1001 char *planMode_env = getenv(
"LAL_FSTAT_FFT_PLAN_MODE" );
1002 char *planGenTimeout_env = getenv(
"LAL_FSTAT_FFT_PLAN_TIMEOUT" );;
1003 int fft_plan_flags = FFTW_MEASURE;
1004 double fft_plan_timeout = FFTW_NO_TIMELIMIT ;
1006 if ( planGenTimeout_env ) {
1008 fft_plan_timeout = strtod( planGenTimeout_env, &
end );
1009 if (
end[0] !=
'\0' ) {
1010 fft_plan_timeout = FFTW_NO_TIMELIMIT;
1014 if ( planMode_env ) {
1015 if ( strcmp( planMode_env,
"ESTIMATE" ) == 0 ) {
1016 fft_plan_flags = FFTW_ESTIMATE;
1019 if ( strcmp( planMode_env,
"MEASURE" ) == 0 ) {
1020 fft_plan_flags = FFTW_MEASURE;
1023 if ( strcmp( planMode_env,
"PATIENT" ) == 0 ) {
1024 fft_plan_flags = FFTW_PATIENT;
1027 *planMode = fft_plan_flags;
1028 *planGenTimeoutSeconds = fft_plan_timeout;
#define GPSSETREAL8(gps, r8)
static int XLALGetFstatTiming_Resamp_intern(const FstatTimingResamp *tiRS, FstatTimingModel *timingModel)
static REAL4 compute_fstat_from_fa_fb(COMPLEX8 Fa, COMPLEX8 Fb, REAL4 A, REAL4 B, REAL4 C, REAL4 E, REAL4 Dinv)
#define LAL_FFTW_WISDOM_UNLOCK
#define LAL_FFTW_WISDOM_LOCK
static int XLALComputeFstatResampGeneric(FstatResults *Fstats, const FstatCommon *common, void *method_data)
static void XLALDestroyResampGenericMethodData(void *method_data)
static int XLALApplySpindownAndFreqShiftGeneric(COMPLEX8 *xOut, const COMPLEX8TimeSeries *xIn, const PulsarDopplerParams *doppler, REAL8 freqShift)
static int XLALComputeFaFb_ResampGeneric(ResampGenericMethodData *resamp, ResampGenericWorkspace *ws, const PulsarDopplerParams thisPoint, REAL8 dFreq, UINT4 numFreqBins, const COMPLEX8TimeSeries *TimeSeries_SRC_a, const COMPLEX8TimeSeries *TimeSeries_SRC_b)
int XLALExtractResampledTimeseries_ResampGeneric(MultiCOMPLEX8TimeSeries **multiTimeSeries_SRC_a, MultiCOMPLEX8TimeSeries **multiTimeSeries_SRC_b, void *method_data)
int XLALGetFstatTiming_ResampGeneric(const void *method_data, FstatTimingGeneric *timingGeneric, FstatTimingModel *timingModel)
static int XLALBarycentricResampleMultiCOMPLEX8TimeSeriesGeneric(ResampGenericMethodData *resamp, const PulsarDopplerParams *thisPoint, const FstatCommon *common)
Performs barycentric resampling on a multi-detector timeseries, updates resampling buffer with result...
int XLALSetupFstatResampGeneric(void **method_data, FstatCommon *common, FstatMethodFuncs *funcs, MultiSFTVector *multiSFTs, const FstatOptionalArgs *optArgs)
static void XLALGetFFTPlanHints(int *planMode, double *planGenTimeoutSeconds)
static void XLALDestroyResampGenericWorkspace(void *workspace)
FstatQuantities
Bit-field of -statistic quantities which can be computed by XLALComputeFstat().
@ FSTATQ_2F
Compute multi-detector .
@ FSTATQ_FAFB_PER_DET
Compute and for each detector.
@ FSTATQ_FAFB
Compute multi-detector and .
@ FSTATQ_NONE
Do not compute -statistic, still compute buffered quantities.
@ FSTATQ_2F_CUDA
Compute multi-detector , store results on CUDA GPU (CUDA implementation of Resamp only).
@ FSTATQ_ATOMS_PER_DET
Compute per-SFT -statistic atoms for each detector (Demod only).
@ FSTATQ_2F_PER_DET
Compute for each detector.
static const REAL8 LAL_FACT_INV[]
void XLALDestroyMultiAMCoeffs(MultiAMCoeffs *multiAMcoef)
Destroy a MultiAMCoeffs structure.
MultiAMCoeffs * XLALComputeMultiAMCoeffs(const MultiDetectorStateSeries *multiDetStates, const MultiNoiseWeights *multiWeights, SkyPosition skypos)
Multi-IFO version of XLALComputeAMCoeffs().
#define XLAL_INIT_DECL(var,...)
void * XLALCalloc(size_t m, size_t n)
void * XLALRealloc(void *p, size_t n)
MultiCOMPLEX8TimeSeries * XLALMultiSFTVectorToCOMPLEX8TimeSeries(const MultiSFTVector *multisfts)
Turn the given multiSFTvector into multiple long COMPLEX8TimeSeries, properly dealing with gaps.
int XLALSincInterpolateCOMPLEX8TimeSeries(COMPLEX8Vector *y_out, const REAL8Vector *t_out, const COMPLEX8TimeSeries *ts_in, UINT4 Dterms)
Interpolate a given regularly-spaced COMPLEX8 timeseries 'ts_in = x_in(j * dt)' onto new samples 'y_o...
void XLALDestroyMultiCOMPLEX8TimeSeries(MultiCOMPLEX8TimeSeries *multiTimes)
Destroy a MultiCOMPLEX8TimeSeries structure.
REAL8 XLALGetCPUTime(void)
#define PULSAR_MAX_DETECTORS
maximal number of detectors we can handle (for static arrays of detector quantities)
#define PULSAR_MAX_SPINS
maximal number of spin-parameters (Freq + spindowns) we can handle
void XLALDestroyMultiSFTVector(MultiSFTVector *multvect)
Destroy a multi SFT-vector.
int XLALMultiSFTVectorResizeBand(MultiSFTVector *multiSFTs, REAL8 f0, REAL8 Band)
Resize the frequency-band of a given multi-SFT vector to [f0, f0+Band].
MultiSSBtimes * XLALGetMultiSSBtimes(const MultiDetectorStateSeries *multiDetStates, SkyPosition skypos, LIGOTimeGPS refTime, SSBprecision precision)
Multi-IFO version of XLALGetSSBtimes().
void XLALDestroyMultiSSBtimes(MultiSSBtimes *multiSSB)
Destroy a MultiSSBtimes structure.
int XLALAddMultiBinaryTimes(MultiSSBtimes **multiSSBOut, const MultiSSBtimes *multiSSBIn, const PulsarDopplerParams *Doppler)
Multi-IFO version of XLALAddBinaryTimes().
int XLALSinCos2PiLUT(REAL4 *sin2pix, REAL4 *cos2pix, REAL8 x)
Calculate sin(2*pi*x) and cos(2*pi*x) to roughly 1e-7 precision using a lookup-table and Taylor-expan...
COORDINATESYSTEM_EQUATORIAL
COMPLEX8TimeSeries * XLALCreateCOMPLEX8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
const LALUnit lalDimensionlessUnit
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
COMPLEX8Vector * XLALCreateCOMPLEX8Vector(UINT4 length)
void XLALDestroyCOMPLEX8Vector(COMPLEX8Vector *vector)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK(assertion,...)
int int XLALPrintWarning(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
This structure contains the per-SFT (weighted) antenna-pattern functions , with the SFT-index,...
REAL4 B
summed antenna-pattern matrix coefficient:
REAL4Vector * b
(weighted) per-SFT antenna-pattern function
REAL4 A
summed antenna-pattern matrix coefficient:
REAL4 C
summed antenna-pattern matrix coefficient:
REAL4Vector * a
(weighted) per-SFT antenna-pattern function
Struct holding the "antenna-pattern" matrix , in terms of the multi-detector scalar product.
REAL4 Dd
determinant factor , such that
MultiDetectorStateSeries * multiDetectorStates
MultiLIGOTimeGPSVector * multiTimestamps
MultiNoiseWeights * multiNoiseWeights
void(* workspace_destroy_func)(void *)
void(* method_data_destroy_func)(void *)
int(* compute_func)(FstatResults *, const FstatCommon *, void *)
Struct of optional 'advanced level' and (potentially method-specific) arguments to be passed to the ...
BOOLEAN resampFFTPowerOf2
Resamp: round up FFT lengths to next power of 2; see FstatMethodType.
BOOLEAN collectTiming
a flag to turn on/off the collection of F-stat-method-specific timing-data
UINT4 Dterms
Number of Dirichlet kernel terms, used by some Demod methods; see FstatMethodType.
XLALComputeFstat() computed results structure.
AntennaPatternMatrix MmunuX[PULSAR_MAX_DETECTORS]
Per detector antenna pattern matrix , used in computing .
AntennaPatternMatrix Mmunu
Antenna pattern matrix , used in computing .
REAL4 * twoFPerDet[PULSAR_MAX_DETECTORS]
If whatWasComputed & FSTATQ_2F_PER_DET is true, the values computed at numFreqBins frequencies space...
FstatQuantities whatWasComputed
Bit-field of which -statistic quantities were computed.
COMPLEX8 * FaPerDet[PULSAR_MAX_DETECTORS]
If whatWasComputed & FSTATQ_PARTS_PER_DET is true, the and values computed at numFreqBins frequenci...
UINT4 numFreqBins
Number of frequencies at which the were computed.
REAL4 * twoF
If whatWasComputed & FSTATQ_2F is true, the multi-detector values computed at numFreqBins frequencie...
COMPLEX8 * Fa
If whatWasComputed & FSTATQ_PARTS is true, the multi-detector and computed at numFreqBins frequenci...
COMPLEX8 * FbPerDet[PULSAR_MAX_DETECTORS]
PulsarDopplerParams doppler
Doppler parameters, including the starting frequency, at which the were computed.
Generic F-stat timing coefficients (times in seconds) [see https://dcc.ligo.org/LIGO-T1600531-v4 for ...
Struct to carry the -statistic method-specific timing model in terms of a list of variable names and...
A vector of 'timestamps' of type LIGOTimeGPS.
REAL8 deltaT
'length' of each timestamp (e.g.
LIGOTimeGPS * data
array of timestamps
UINT4 length
number of timestamps
Multi-IFO container for antenna-pattern coefficients and atenna-pattern matrix .
AMCoeffs ** data
noise-weighted AM-coeffs , and
AntennaPatternMatrix Mmunu
antenna-pattern matrix
Multi-IFO container for COMPLEX8 resampled timeseries.
COMPLEX8TimeSeries ** data
array of COMPLEX8TimeSeries (pointers)
UINT4 length
number of IFOs
LIGOTimeGPSVector ** data
timestamps vector for each ifo
A collection of SFT vectors – one for each IFO in a multi-IFO search.
Multi-IFO container for SSB timings.
SSBtimes ** data
array of SSBtimes (pointers)
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
REAL8 period
Binary: orbital period (sec)
LIGOTimeGPS tp
Binary: time of observed periapsis passage (in SSB)
PulsarSpins fkdot
Intrinsic spins: [Freq, f1dot, f2dot, ...
REAL8 Delta
Sky position: DEC (latitude) in equatorial coords and radians.
LIGOTimeGPS refTime
Reference time of pulsar parameters (in SSB!)
REAL8 Alpha
Sky position: RA (longitude) in equatorial coords and radians.
REAL8 ecc
Binary: orbital eccentricity.
REAL8 asini
Binary: projected, normalized orbital semi-major axis (s).
REAL8 argp
Binary: argument of periapsis (radians)
AntennaPatternMatrix Mmunu
MultiSSBtimes * multiSSBtimes
FstatTimingResamp timingResamp
MultiSSBtimes * multiBinaryTimes
MultiCOMPLEX8TimeSeries * multiTimeSeries_DET
MultiAMCoeffs * multiAMcoef
PulsarDopplerParams prev_doppler
MultiCOMPLEX8TimeSeries * multiTimeSeries_SRC_a
AntennaPatternMatrix MmunuX[PULSAR_MAX_DETECTORS]
MultiCOMPLEX8TimeSeries * multiTimeSeries_SRC_b
FstatTimingGeneric timingGeneric
REAL8Vector * SRCtimes_DET
COMPLEX8Vector * TStmp1_SRC
COMPLEX8Vector * TStmp2_SRC
Simple container for two REAL8-vectors, namely the SSB-timings DeltaT_alpha and Tdot_alpha,...
REAL8Vector * Tdot
dT/dt : time-derivative of SSB-time wrt local time for SFT-alpha
REAL8Vector * DeltaT
Time-difference of SFT-alpha - tau0 in SSB-frame.
LIGOTimeGPS refTime
reference-time 'tau0'