Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
ComputeFstatBenchmark.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2017 Maximillian Bensch
3* Copyright (C) 2015 Reinhard Prix
4*
5* This program is free software; you can redistribute it and/or modify
6* it under the terms of the GNU General Public License as published by
7* the Free Software Foundation; either version 2 of the License, or
8* (at your option) any later version.
9*
10* This program is distributed in the hope that it will be useful,
11* but WITHOUT ANY WARRANTY; without even the implied warranty of
12* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13* GNU General Public License for more details.
14*
15* You should have received a copy of the GNU General Public License
16* along with with program; see the file COPYING. If not, write to the
17* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
18* MA 02110-1301 USA
19*/
20
21#include "config.h"
22
23#include <lal/XLALError.h>
24#include <lal/LALBarycenter.h>
25#include <lal/LALInitBarycenter.h>
26#include <lal/LogPrintf.h>
27#include <lal/CWMakeFakeData.h>
28#include <lal/LALConstants.h>
29#include <lal/ExtrapolatePulsarSpins.h>
30#include <lal/ComputeFstat.h>
31#include <lal/DetectorStates.h>
32#include <lal/LFTandTSutils.h>
33#include <lal/LALString.h>
34#include <lal/UserInput.h>
35#include <lal/LALPulsarVCSInfo.h>
36
37// benchmark ComputeFstat() functions for performance and memory usage
39
40typedef struct {
41 int FstatMethod; //!< select which method/algorithm to use to compute the F-statistic
44 REAL8Range Freq;
46 REAL8Range f2dot;
47 REAL8Range FreqResolution;
48 REAL8Range orbitasini;
49 REAL8Range orbitPeriod;
50 REAL8Range orbitEcc;
51 REAL8Range orbitArgp;
52 LIGOTimeGPSRange orbitTp;
53 INT4Range numFreqBins;
54 INT4Range Tseg;
55 INT4 numSegments;
57 CHAR *outputInfo;
58 INT4 numTrials;
60
61 // ----- developer options
62 CHAR *ephemEarth; /**< Earth ephemeris file to use */
63 CHAR *ephemSun; /**< Sun ephemeris file to use */
64
65 REAL8 Tsft;
66 BOOLEAN sharedWorkspace; // useful for checking workspace sharing for Resampling
67 BOOLEAN perSegmentSFTs; // Weave vs GCT convention: GCT loads SFT frequency ranges globally, Weave loads them per segment (more efficient)
68 BOOLEAN resampFFTPowerOf2;
69 INT4 Dterms;
71
72 BOOLEAN version; // output code version
74
75// ---------- main ----------
76int
77main( int argc, char *argv[] )
78{
79
80 CHAR *VCSInfoString;
81 XLAL_CHECK_MAIN( ( VCSInfoString = XLALVCSInfoString( lalPulsarVCSInfoList, 0, "%% " ) ) != NULL, XLAL_EFUNC );
82
83 // ---------- handle user input ----------
85 UserInput_t *uvar = &uvar_s;
86
87 uvar->FstatMethod = FMETHOD_RESAMP_BEST;
88 uvar->Alpha[0] = 0;
89 uvar->Alpha[1] = LAL_TWOPI;
90 uvar->Delta[0] = -LAL_PI / 2;
91 uvar->Delta[1] = LAL_PI / 2;
92 uvar->Freq[0] = 100;
93 uvar->Freq[1] = 1000;
94 uvar->f1dot[0] = -3e-9;
95 uvar->f1dot[1] = 0;
96 uvar->f2dot[0] = 0;
97 uvar->f2dot[1] = 1e-16;
98 uvar->FreqResolution[0] = 0.1;
99 uvar->FreqResolution[1] = 1;
100 uvar->numFreqBins[0] = 1000;
101 uvar->numFreqBins[1] = 100000;
102 uvar->Tseg[0] = 10 * 3600;
103 uvar->Tseg[1] = 250 * 3600;
104 uvar->orbitasini[0] = 0;
105 uvar->orbitasini[1] = 0;
106 uvar->orbitPeriod[0] = 0;
107 uvar->orbitPeriod[1] = 0;
108 uvar->orbitEcc[0] = 0;
109 uvar->orbitEcc[1] = 0;
110 uvar->orbitArgp[0] = 0;
111 uvar->orbitArgp[1] = 0;
112 uvar->orbitTp[0].gpsSeconds = 0;
113 uvar->orbitTp[1].gpsSeconds = 0;
114
115 uvar->numSegments = 90;
116 uvar->numTrials = 1;
117 uvar->startTime.gpsSeconds = 711595934;
118 uvar->Tsft = 1800;
119 uvar->sharedWorkspace = 1;
120 uvar->resampFFTPowerOf2 = FstatOptionalArgsDefaults.resampFFTPowerOf2;
121 uvar->perSegmentSFTs = 1;
122
123 uvar->Dterms = FstatOptionalArgsDefaults.Dterms;
124
125 uvar->ephemEarth = XLALStringDuplicate( "earth00-40-DE405.dat.gz" );
126 uvar->ephemSun = XLALStringDuplicate( "sun00-40-DE405.dat.gz" );
127 uvar->randSeed = 1;
128
129 XLAL_CHECK_MAIN( ( uvar->IFOs = XLALCreateStringVector( "H1", NULL ) ) != NULL, XLAL_EFUNC );
130 uvar->outputInfo = NULL;
131
132 XLAL_CHECK( XLALRegisterUvarAuxDataMember( FstatMethod, UserEnum, XLALFstatMethodChoices(), 0, OPTIONAL, "F-statistic method to use" ) == XLAL_SUCCESS, XLAL_EFUNC );
133 XLAL_CHECK_MAIN( XLALRegisterUvarMember( Alpha, RAJRange, 0, OPTIONAL, "Skyposition [drawn isotropically]: Range in 'Alpha' = right ascension)" ) == XLAL_SUCCESS, XLAL_EFUNC );
134 XLAL_CHECK_MAIN( XLALRegisterUvarMember( Delta, DECJRange, 0, OPTIONAL, "Skyposition [drawn isotropically]: Range in 'Delta' = declination" ) == XLAL_SUCCESS, XLAL_EFUNC );
135 XLAL_CHECK_MAIN( XLALRegisterUvarMember( Freq, REAL8Range, 0, OPTIONAL, "Search frequency in Hz [range to draw from]" ) == XLAL_SUCCESS, XLAL_EFUNC );
136 XLAL_CHECK_MAIN( XLALRegisterUvarMember( f1dot, REAL8Range, 0, OPTIONAL, "Search 1st spindown f1dot in Hz/s [range to draw from]" ) == XLAL_SUCCESS, XLAL_EFUNC );
137 XLAL_CHECK_MAIN( XLALRegisterUvarMember( f2dot, REAL8Range, 0, OPTIONAL, "Search 2nd spindown f2dot in Hz/s^2 [range to draw from]" ) == XLAL_SUCCESS, XLAL_EFUNC );
138 XLAL_CHECK_MAIN( XLALRegisterUvarMember( FreqResolution, REAL8Range, 0, OPTIONAL, "Frequency resolution 'R' in natural units 1/Tseg such that: dFreq = R/Tseg) [range to draw from]" ) == XLAL_SUCCESS, XLAL_EFUNC );
139
140 XLAL_CHECK_MAIN( XLALRegisterUvarMember( orbitasini, REAL8Range, 0, OPTIONAL, "Binary Orbit: Projected semi-major axis in light-seconds [range to draw from]" ) == XLAL_SUCCESS, XLAL_EFUNC );
141 XLAL_CHECK_MAIN( XLALRegisterUvarMember( orbitPeriod, REAL8Range, 0, OPTIONAL, "Binary Orbit: Period in seconds [range to draw from]" ) == XLAL_SUCCESS, XLAL_EFUNC );
142 XLAL_CHECK_MAIN( XLALRegisterUvarMember( orbitTp, EPOCHRange, 0, OPTIONAL, "Binary Orbit: (true) epoch of periapsis: use 'xx.yy[GPS|MJD]' format. [range to draw from]" ) == XLAL_SUCCESS, XLAL_EFUNC );
143 XLAL_CHECK_MAIN( XLALRegisterUvarMember( orbitArgp, REAL8Range, 0, OPTIONAL, "Binary Orbit: Orbital argument of periapse in radians [range to draw from]" ) == XLAL_SUCCESS, XLAL_EFUNC );
144 XLAL_CHECK_MAIN( XLALRegisterUvarMember( orbitEcc, REAL8Range, 0, OPTIONAL, "Binary Orbit: Orbital eccentricity [range to draw from]" ) == XLAL_SUCCESS, XLAL_EFUNC );
145
146 XLAL_CHECK_MAIN( XLALRegisterUvarMember( startTime, EPOCH, 0, OPTIONAL, "Start time of first segment" ) == XLAL_SUCCESS, XLAL_EFUNC );
147 XLAL_CHECK_MAIN( XLALRegisterUvarMember( Tseg, INT4Range, 0, OPTIONAL, "Coherent segment length in seconds [range to draw from]" ) == XLAL_SUCCESS, XLAL_EFUNC );
148 XLAL_CHECK_MAIN( XLALRegisterUvarMember( numSegments, INT4, 0, OPTIONAL, "Number of semi-coherent segments" ) == XLAL_SUCCESS, XLAL_EFUNC );
149 XLAL_CHECK_MAIN( XLALRegisterUvarMember( numFreqBins, INT4Range, 0, OPTIONAL, "Number of frequency bins to search [range to draw from]" ) == XLAL_SUCCESS, XLAL_EFUNC );
150 XLAL_CHECK_MAIN( XLALRegisterUvarMember( IFOs, STRINGVector, 0, OPTIONAL, "IFOs to use [list]" ) == XLAL_SUCCESS, XLAL_EFUNC );
151 XLAL_CHECK_MAIN( XLALRegisterUvarMember( numTrials, INT4, 0, OPTIONAL, "Number of repeated trials to run (with randomized parameters drawn from ranges)" ) == XLAL_SUCCESS, XLAL_EFUNC );
152
153 XLAL_CHECK_MAIN( XLALRegisterUvarMember( sharedWorkspace, BOOLEAN, 0, OPTIONAL, "Use workspace sharing across segments (only used in Resampling)" ) == XLAL_SUCCESS, XLAL_EFUNC );
154 XLAL_CHECK_MAIN( XLALRegisterUvarMember( perSegmentSFTs, BOOLEAN, 0, OPTIONAL, "Weave vs GCT: GCT determines and loads SFT frequency ranges globally, Weave does that per segment (more efficient)" ) == XLAL_SUCCESS, XLAL_EFUNC );
155 XLAL_CHECK_MAIN( XLALRegisterUvarMember( resampFFTPowerOf2, BOOLEAN, 0, OPTIONAL, "For Resampling methods: enforce FFT length to be a power of two (by rounding up)" ) == XLAL_SUCCESS, XLAL_EFUNC );
156
157 XLAL_CHECK_MAIN( XLALRegisterUvarMember( Dterms, INT4, 0, OPTIONAL, "Number of kernel terms (single-sided) in\na) Dirichlet kernel if FstatMethod=Demod*\nb) sinc-interpolation if FstatMethod=Resamp*" ) == XLAL_SUCCESS, XLAL_EFUNC );
158
159 XLAL_CHECK_MAIN( XLALRegisterUvarMember( outputInfo, STRING, 0, OPTIONAL, "Append Resampling internal info into this file" ) == XLAL_SUCCESS, XLAL_EFUNC );
160
161 XLAL_CHECK_MAIN( XLALRegisterUvarMember( Tsft, REAL8, 0, DEVELOPER, "SFT length" ) == XLAL_SUCCESS, XLAL_EFUNC );
162 XLAL_CHECK_MAIN( XLALRegisterUvarMember( ephemEarth, STRING, 0, DEVELOPER, "Earth ephemeris file to use" ) == XLAL_SUCCESS, XLAL_EFUNC );
163 XLAL_CHECK_MAIN( XLALRegisterUvarMember( ephemSun, STRING, 0, DEVELOPER, "Sun ephemeris file to use" ) == XLAL_SUCCESS, XLAL_EFUNC );
164 XLAL_CHECK_MAIN( XLALRegisterUvarMember( randSeed, INT4, 0, DEVELOPER, "Random seed to use for rand() to draw randomized parameters." ) == XLAL_SUCCESS, XLAL_EFUNC );
165
166 BOOLEAN should_exit = 0;
168 if ( should_exit ) {
169 return EXIT_FAILURE;
170 }
171
172 BOOLEAN have_orbitasini = XLALUserVarWasSet( &uvar->orbitasini );
173 BOOLEAN have_orbitPeriod = XLALUserVarWasSet( &uvar->orbitPeriod );
174 BOOLEAN have_orbitTp = XLALUserVarWasSet( &uvar->orbitTp );
175 BOOLEAN have_orbitEcc = XLALUserVarWasSet( &uvar->orbitEcc );
176 BOOLEAN have_orbitArgp = XLALUserVarWasSet( &uvar->orbitArgp );
177
178 if ( have_orbitasini || have_orbitEcc || have_orbitPeriod || have_orbitArgp || have_orbitTp ) {
179 XLAL_CHECK( uvar->orbitasini[0] >= 0, XLAL_EDOM );
180 XLAL_CHECK( ( uvar->orbitasini[1] == 0 ) || ( have_orbitPeriod && ( uvar->orbitPeriod[0] > 0 ) && have_orbitTp ), XLAL_EINVAL, "If orbitasini>0 then we also need 'orbitPeriod>0' and 'orbitTp'\n" );
181 XLAL_CHECK( ( uvar->orbitEcc[0] >= 0 ) && ( uvar->orbitEcc[1] <= 1 ), XLAL_EDOM );
182 }
183
184 // produce log-string (for output-file headers)
185 CHAR *logstring = NULL;
186 XLAL_CHECK_MAIN( ( logstring = XLALStringAppend( logstring, "%% cmdline: " ) ) != NULL, XLAL_EFUNC );
187 CHAR *cmdline;
189 XLAL_CHECK_MAIN( ( logstring = XLALStringAppend( logstring, cmdline ) ) != NULL, XLAL_EFUNC );
190 XLALFree( cmdline );
191 XLAL_CHECK_MAIN( ( logstring = XLALStringAppend( logstring, "\n" ) ) != NULL, XLAL_EFUNC );
192 XLAL_CHECK_MAIN( ( logstring = XLALStringAppend( logstring, VCSInfoString ) ) != NULL, XLAL_EFUNC );
193
194 // check user input
195 XLAL_CHECK_MAIN( uvar->numSegments >= 1, XLAL_EINVAL );
196 XLAL_CHECK_MAIN( uvar->Tsft > 1, XLAL_EINVAL );
197 XLAL_CHECK_MAIN( uvar->numTrials >= 1, XLAL_EINVAL );
198 // ---------- end: handle user input ----------
199 srand( uvar->randSeed ); // set random seed
200
201 // common setup over repeated trials
202 EphemerisData *ephem;
203 XLAL_CHECK_MAIN( ( ephem = XLALInitBarycenter( uvar->ephemEarth, uvar->ephemSun ) ) != NULL, XLAL_EFUNC );
205
206 UINT4 numDetectors = uvar->IFOs->length;
207 // ----- setup injection and data parameters
208 LIGOTimeGPSVector *startTime_l, *endTime_l;
209 XLAL_CHECK_MAIN( ( startTime_l = XLALCreateTimestampVector( uvar->numSegments ) ) != NULL, XLAL_EFUNC );
210 XLAL_CHECK_MAIN( ( endTime_l = XLALCreateTimestampVector( uvar->numSegments ) ) != NULL, XLAL_EFUNC );
211
212 // ----- setup optional Fstat arguments
214 MultiNoiseFloor XLAL_INIT_DECL( injectSqrtSX );
215 injectSqrtSX.length = numDetectors;
216 for ( UINT4 X = 0; X < numDetectors; X ++ ) {
217 injectSqrtSX.sqrtSn[X] = 1;
218 }
219 optionalArgs.injectSqrtSX = &injectSqrtSX;
220 optionalArgs.FstatMethod = uvar->FstatMethod;
221 optionalArgs.collectTiming = 1;
222 optionalArgs.resampFFTPowerOf2 = uvar->resampFFTPowerOf2;
223 optionalArgs.Dterms = uvar->Dterms;
224
225 FILE *timingLogFILE = NULL;
226 FILE *timingParFILE = NULL;
227 if ( uvar->outputInfo != NULL ) {
228 char *parFname = NULL;
229 XLAL_CHECK_MAIN( ( parFname = XLALStringAppend( parFname, uvar->outputInfo ) ) != NULL, XLAL_EFUNC );
230 XLAL_CHECK_MAIN( ( parFname = XLALStringAppend( parFname, ".pars" ) ) != NULL, XLAL_EFUNC );
231
232 XLAL_CHECK_MAIN( ( timingLogFILE = fopen( uvar->outputInfo, "ab" ) ) != NULL, XLAL_ESYS, "Failed to open '%s' for appending\n", uvar->outputInfo );
233 XLAL_CHECK_MAIN( ( timingParFILE = fopen( parFname, "ab" ) ) != NULL, XLAL_ESYS, "Failed to open '%s' for appending\n", parFname );
234 XLALFree( parFname );
235 fprintf( timingLogFILE, "%s\n", logstring );
236 fprintf( timingParFILE, "%s\n", logstring );
237 fprintf( timingParFILE, "%%%%%8s %20s %20s %20s %20s %20s %20s %20s %20s %12s %20s %20s %20s %20s %20s\n",
238 "Nseg", "Tseg", "Freq", "FreqBand", "dFreq", "f1dot", "f2dot", "Alpha", "Delta", "memUsageMB", "asini", "period", "ecc", "argp", "tp" );
239 }
240 FstatInputVector *inputs;
241 FstatQuantities whatToCompute = ( FSTATQ_2F | FSTATQ_2F_PER_DET );
242 FstatResults *results = NULL;
243
244#define drawFromREAL8Range(range) (range[0] + (range[1] - range[0]) * rand() / RAND_MAX )
245#define drawFromINT4Range(range) (range[0] + (INT4)round(1.0*(range[1] - range[0]) * rand() / RAND_MAX) )
246#define drawFromEPOCHRange(epoch, range) XLALGPSSetREAL8 ( epoch, (XLALGPSGetREAL8(&range[0]) + round(1.0*(XLALGPSGetREAL8(&range[1]) - XLALGPSGetREAL8(&range[0])) * rand() / RAND_MAX) ))
247 // ---------- main loop over repeated trials: randomize uniformly over input ranges ----------
248 for ( INT4 i = 0; i < uvar->numTrials; i ++ ) {
249 UINT4 Tseg_i = drawFromINT4Range( uvar->Tseg );
250 Tseg_i = ( UINT4 ) uvar->Tsft * ceil( Tseg_i / uvar->Tsft );
251 SFTCatalog **catalogs;
252 XLAL_CHECK_MAIN( ( catalogs = XLALCalloc( uvar->numSegments, sizeof( catalogs[0] ) ) ) != NULL, XLAL_ENOMEM );
253 for ( INT4 l = 0; l < uvar->numSegments; l ++ ) {
254 startTime_l->data[l] = ( l == 0 ) ? uvar->startTime : endTime_l->data[l - 1];
255 endTime_l->data[l] = startTime_l->data[l];
256 endTime_l->data[l].gpsSeconds += Tseg_i;
257
259 XLAL_CHECK_MAIN( ( multiTimestamps = XLALMakeMultiTimestamps( startTime_l->data[l], Tseg_i, uvar->Tsft, 0, numDetectors ) ) != NULL, XLAL_EFUNC );
260 XLAL_CHECK_MAIN( ( catalogs[l] = XLALMultiAddToFakeSFTCatalog( NULL, uvar->IFOs, multiTimestamps ) ) != NULL, XLAL_EFUNC );
262 } // for l < numSegments
263
264 PulsarSpinRange XLAL_INIT_DECL( spinRange_i );
265
266 LIGOTimeGPS refTime = { uvar->startTime.gpsSeconds + 0.5 * uvar->numSegments * Tseg_i, 0 };
267 spinRange_i.refTime = refTime;
268 spinRange_i.fkdot[0] = drawFromREAL8Range( uvar->Freq );
269 spinRange_i.fkdot[1] = drawFromREAL8Range( uvar->f1dot );
270 spinRange_i.fkdot[2] = drawFromREAL8Range( uvar->f2dot );
271
273 Doppler_i.refTime = refTime;
274 Doppler_i.Alpha = drawFromREAL8Range( uvar->Alpha );
275 REAL8Range sDeltaRange;
276 sDeltaRange[0] = sin( uvar->Delta[0] );
277 sDeltaRange[1] = sin( uvar->Delta[1] );
278 Doppler_i.Delta = asin( drawFromREAL8Range( sDeltaRange ) );
279 memcpy( &Doppler_i.fkdot, &spinRange_i.fkdot, sizeof( Doppler_i.fkdot ) );;
280 Doppler_i.period = drawFromREAL8Range( uvar->orbitPeriod );
281 Doppler_i.ecc = drawFromREAL8Range( uvar->orbitEcc );
282 Doppler_i.asini = drawFromREAL8Range( uvar->orbitasini );
283 drawFromEPOCHRange( &Doppler_i.tp, uvar->orbitTp );
284 Doppler_i.argp = drawFromREAL8Range( uvar->orbitArgp );
285
286 UINT4 numFreqBins_i = drawFromINT4Range( uvar->numFreqBins );
287 REAL8 FreqResolution_i = drawFromREAL8Range( uvar->FreqResolution );
288
289 REAL8 dFreq_i = FreqResolution_i / Tseg_i;
290 REAL8 FreqBand_i = numFreqBins_i * dFreq_i;
291
292 XLAL_CHECK_MAIN( ( inputs = XLALCreateFstatInputVector( uvar->numSegments ) ) != NULL, XLAL_EFUNC );
293
294 fprintf( stderr, "trial %d/%d: Tseg = %.1f d, numSegments = %d, Alpha = %.2f rad, Delta = %.2f rad, Freq = %.6f Hz, f1dot = %.1e Hz/s, f2dot = %.1e Hz/s^2, R = %.2f, numFreqBins = %d, asini = %.2f, period = %.2f, ecc = %.2f, argp = %.2f, tp=%"LAL_GPS_FORMAT" [dFreq = %.2e Hz, FreqBand = %.2e Hz]\n",
295 i + 1, uvar->numTrials, Tseg_i / 86400.0, uvar->numSegments, Doppler_i.Alpha, Doppler_i.Delta, Doppler_i.fkdot[0], Doppler_i.fkdot[1], Doppler_i.fkdot[2], FreqResolution_i, numFreqBins_i, Doppler_i.asini, Doppler_i.period, Doppler_i.ecc, Doppler_i.argp, LAL_GPS_PRINT( Doppler_i.tp ), dFreq_i, FreqBand_i );
296
297 spinRange_i.fkdotBand[0] = FreqBand_i;
298 REAL8 minCoverFreq_il, maxCoverFreq_il;
299 // GCT convention: determine global SFT frequency band for all segments
300 if ( ! uvar->perSegmentSFTs ) {
301 XLAL_CHECK_MAIN( XLALCWSignalCoveringBand( &minCoverFreq_il, &maxCoverFreq_il, &startTime_l->data[0], &endTime_l->data[uvar->numSegments - 1], &spinRange_i, Doppler_i.asini, Doppler_i.period, Doppler_i.ecc ) == XLAL_SUCCESS, XLAL_EFUNC );
302 }
303 // create per-segment input structs
304 for ( INT4 l = 0; l < uvar->numSegments; l ++ ) {
305 if ( uvar->sharedWorkspace && l > 0 ) {
306 optionalArgs.prevInput = inputs->data[0];
307 } else {
308 optionalArgs.prevInput = NULL;
309 }
310 // Weave convention: determine per-segment SFT frequency band
311 if ( uvar->perSegmentSFTs ) {
312 XLAL_CHECK_MAIN( XLALCWSignalCoveringBand( &minCoverFreq_il, &maxCoverFreq_il, &startTime_l->data[l], &endTime_l->data[l], &spinRange_i, Doppler_i.asini, Doppler_i.period, Doppler_i.ecc ) == XLAL_SUCCESS, XLAL_EFUNC );
313 // fprintf ( stderr, "minCoverFreq = %.16g, maxCoverFreq = %.16g: startTime = %d, Tseg = %d, FreqMax = %.16g, f1dot = %.16g, f2dot = %.16g, refTime = %d = %d + %d\n",
314 // minCoverFreq_il, maxCoverFreq_il, startTime_l->data[l].gpsSeconds, Tseg_i, spinRange_i.fkdot[0] + spinRange_i.fkdotBand[0], spinRange_i.fkdot[1], spinRange_i.fkdot[2],
315 // refTime.gpsSeconds, startTime_l->data[l].gpsSeconds, refTime.gpsSeconds - startTime_l->data[l].gpsSeconds );
316 }
317 XLAL_CHECK_MAIN( ( inputs->data[l] = XLALCreateFstatInput( catalogs[l], minCoverFreq_il, maxCoverFreq_il, dFreq_i, ephem, &optionalArgs ) ) != NULL, XLAL_EFUNC );
318 }
319 for ( INT4 l = 0; l < uvar->numSegments; l ++ ) {
320 XLALDestroySFTCatalog( catalogs[l] );
321 }
322 XLALFree( catalogs );
323
324 // ----- compute Fstatistics over segments
325 for ( INT4 l = 0; l < uvar->numSegments; l ++ ) {
326 XLAL_CHECK_MAIN( XLALComputeFstat( &results, inputs->data[l], &Doppler_i, numFreqBins_i, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
327
328 // ----- output timing details to file if requested
329 if ( timingLogFILE != NULL ) {
330 XLAL_CHECK_MAIN( XLALAppendFstatTiming2File( inputs->data[l], timingLogFILE, ( l == 0 ) && ( i == 0 ) ) == XLAL_SUCCESS, XLAL_EFUNC );
331 }
332 } // for l < numSegments
333
335 REAL8 memUsage = memEnd - memBase;
336 const char *FmethodName = XLALGetFstatInputMethodName( inputs->data[0] );
337 fprintf( stderr, "%-15s: memoryUsage = %6.1f MB\n", FmethodName, memUsage );
338
339 if ( timingParFILE != NULL ) {
340 fprintf( timingParFILE, "%10d %20d %20.16g %20.16g %20.16g %20.16g %20.16g %20.16g %20.16g %12g %20.16g %20.16g %20.16g %20.16g %"LAL_GPS_FORMAT"\n",
341 uvar->numSegments, Tseg_i, Doppler_i.fkdot[0], FreqBand_i, dFreq_i, Doppler_i.fkdot[1], Doppler_i.fkdot[2], Doppler_i.Alpha, Doppler_i.Delta, memUsage, Doppler_i.asini, Doppler_i.period, Doppler_i.ecc, Doppler_i.argp, LAL_GPS_PRINT( Doppler_i.tp )
342 );
343 }
344
346 } // for i < numTrials
347
348 // ----- free memory ----------
349 if ( timingLogFILE != NULL ) {
350 fclose( timingLogFILE );
351 }
352 if ( timingParFILE != NULL ) {
353 fclose( timingParFILE );
354 }
355
356 XLALDestroyFstatResults( results );
359 XLALFree( VCSInfoString );
360 XLALFree( logstring );
361 XLALDestroyTimestampVector( startTime_l );
362 XLALDestroyTimestampVector( endTime_l );
363
365
366 return XLAL_SUCCESS;
367
368} // main()
369
370
371// --------------------------------------------------------------------------------
372// code to read current process RSS memory usage from /proc, taken from
373// https://stackoverflow.com/questions/63166/how-to-determine-cpu-and-memory-consumption-from-inside-a-process
374static int parseLine( char *line );
375
376REAL8
378{
379 FILE *file = fopen( "/proc/self/status", "r" );
380 int result = -1;
381 char line[128];
382 if ( file == NULL ) {
383 return result;
384 }
385 while ( fgets( line, sizeof( line ), file ) != NULL ) {
386 if ( strncmp( line, "VmRSS:", 6 ) == 0 ) {
387 result = parseLine( line );
388 break;
389 }
390 }
391 fclose( file );
392 return ( result / 1024.0 );
393} // XLALGetCurrentHeapUsageMB()
394
395static int parseLine( char *line )
396{
397 // This assumes that a digit will be found and the line ends in " Kb".
398 int i = strlen( line );
399 const char *p = line;
400 while ( *p < '0' || *p > '9' ) {
401 p++;
402 }
403 line[i - 3] = '\0';
404 i = atoi( p );
405 return i;
406}
407// --------------------------------------------------------------------------------
int main(int argc, char *argv[])
REAL8 XLALGetCurrentHeapUsageMB(void)
static int parseLine(char *line)
#define drawFromEPOCHRange(epoch, range)
#define drawFromINT4Range(range)
#define drawFromREAL8Range(range)
void LALCheckMemoryLeaks(void)
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
#define STRING(a)
#define fprintf
int l
double e
const CHAR * XLALGetFstatInputMethodName(const FstatInput *input)
Returns the human-readable name of the -statistic method being used by a FstatInput structure.
Definition: ComputeFstat.c:672
const FstatOptionalArgs FstatOptionalArgsDefaults
Global initializer for setting FstatOptionalArgs to default values.
Definition: ComputeFstat.c:90
const UserChoices * XLALFstatMethodChoices(void)
Return pointer to a static array of all (available) FstatMethodType choices.
void XLALDestroyFstatInputVector(FstatInputVector *inputs)
Free all memory associated with a FstatInputVector structure.
Definition: ComputeFstat.c:252
int XLALAppendFstatTiming2File(const FstatInput *input, FILE *fp, BOOLEAN printHeader)
FstatInputVector * XLALCreateFstatInputVector(const UINT4 length)
Create a FstatInputVector of the given length, for example for setting up F-stat searches over severa...
Definition: ComputeFstat.c:231
FstatQuantities
Bit-field of -statistic quantities which can be computed by XLALComputeFstat().
Definition: ComputeFstat.h:94
FstatInput * XLALCreateFstatInput(const SFTCatalog *SFTcatalog, const REAL8 minCoverFreq, const REAL8 maxCoverFreq, const REAL8 dFreq, const EphemerisData *ephemerides, const FstatOptionalArgs *optionalArgs)
Create a fully-setup FstatInput structure for computing the -statistic using XLALComputeFstat().
Definition: ComputeFstat.c:359
int XLALComputeFstat(FstatResults **Fstats, FstatInput *input, const PulsarDopplerParams *doppler, const UINT4 numFreqBins, const FstatQuantities whatToCompute)
Compute the -statistic over a band of frequencies.
Definition: ComputeFstat.c:760
void XLALDestroyFstatResults(FstatResults *Fstats)
Free all memory associated with a FstatResults structure.
Definition: ComputeFstat.c:934
@ FMETHOD_RESAMP_BEST
Resamp: best guess of the fastest available implementation
Definition: ComputeFstat.h:125
@ FSTATQ_2F
Compute multi-detector .
Definition: ComputeFstat.h:96
@ FSTATQ_2F_PER_DET
Compute for each detector.
Definition: ComputeFstat.h:98
#define LAL_GPS_FORMAT
#define LAL_GPS_PRINT(gps)
int XLALCWSignalCoveringBand(REAL8 *minCoverFreq, REAL8 *maxCoverFreq, const LIGOTimeGPS *time1, const LIGOTimeGPS *time2, const PulsarSpinRange *spinRange, const REAL8 binaryMaxAsini, const REAL8 binaryMinPeriod, const REAL8 binaryMaxEcc)
Determines a frequency band which covers the frequency evolution of a band of CW signals between two ...
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
#define LAL_PI
#define LAL_TWOPI
unsigned char BOOLEAN
double REAL8
#define XLAL_INIT_DECL(var,...)
char CHAR
uint32_t UINT4
int32_t INT4
void * XLALCalloc(size_t m, size_t n)
void XLALFree(void *p)
char char * XLALStringDuplicate(const char *s)
int char * XLALStringAppend(char *s, const char *append)
char * XLALVCSInfoString(const LALVCSInfoList vcs_list, const int verbose, const char *prefix)
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
Definition: SFTcatalog.c:329
MultiLIGOTimeGPSVector * XLALMakeMultiTimestamps(LIGOTimeGPS tStart, REAL8 Tspan, REAL8 Tsft, REAL8 Toverlap, UINT4 numDet)
Same as XLALMakeTimestamps() just for several detectors, additionally specify the number of detectors...
SFTCatalog * XLALMultiAddToFakeSFTCatalog(SFTCatalog *catalog, const LALStringVector *detectors, const MultiLIGOTimeGPSVector *timestamps)
Multi-detector and multi-timestamp wrapper of XLALAddToFakeSFTCatalog().
Definition: SFTcatalog.c:749
LIGOTimeGPSVector * XLALCreateTimestampVector(UINT4 len)
Allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:47
void XLALDestroyMultiTimestamps(MultiLIGOTimeGPSVector *multiTS)
Destroy a MultiLIGOTimeGPSVector timestamps vector.
void XLALDestroyTimestampVector(LIGOTimeGPSVector *vect)
De-allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:69
LALStringVector * XLALCreateStringVector(const CHAR *str1,...)
int XLALUserVarReadAllInput(BOOLEAN *should_exit, int argc, char *argv[], const LALVCSInfoList vcs_list)
void XLALDestroyUserVars(void)
#define XLALRegisterUvarMember(name, type, option, category,...)
void CHAR * XLALUserVarGetLog(UserVarLogFormat format)
int XLALUserVarWasSet(const void *cvar)
#define XLALRegisterUvarAuxDataMember(name, type, cdata, option, category,...)
UVAR_LOGFMT_CMDLINE
LIGOTimeGPS LIGOTimeGPSRange[2]
REAL8 REAL8Range[2]
INT4 INT4Range[2]
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_MAIN(assertion,...)
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EDOM
XLAL_ESYS
XLAL_EINVAL
size_t numDetectors
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
A vector of XLALComputeFstat() input data structures, for e.g.
Definition: ComputeFstat.h:82
FstatInput ** data
Pointer to the data array.
Definition: ComputeFstat.h:87
Struct of optional 'advanced level' and (potentially method-specific) arguments to be passed to the ...
Definition: ComputeFstat.h:137
MultiNoiseFloor * injectSqrtSX
Single-sided PSD values for fake Gaussian noise to be added to SFT data.
Definition: ComputeFstat.h:144
BOOLEAN resampFFTPowerOf2
Resamp: round up FFT lengths to next power of 2; see FstatMethodType.
Definition: ComputeFstat.h:148
FstatInput * prevInput
An FstatInput structure from a previous call to XLALCreateFstatInput(); may contain common workspace ...
Definition: ComputeFstat.h:146
BOOLEAN collectTiming
a flag to turn on/off the collection of F-stat-method-specific timing-data
Definition: ComputeFstat.h:147
UINT4 Dterms
Number of Dirichlet kernel terms, used by some Demod methods; see FstatMethodType.
Definition: ComputeFstat.h:140
FstatMethodType FstatMethod
Method to use for computing the -statistic.
Definition: ComputeFstat.h:142
XLALComputeFstat() computed results structure.
Definition: ComputeFstat.h:202
A vector of 'timestamps' of type LIGOTimeGPS.
Definition: SFTfileIO.h:188
LIGOTimeGPS * data
array of timestamps
Definition: SFTfileIO.h:193
A collection of (multi-IFO) LIGOTimeGPSVector time-stamps vectors.
Definition: SFTfileIO.h:198
array of detector-specific 'noise floors' (ie PSD values), assumed constant over the frequency-band o...
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
Contains a "spin-range", ie spins and corresponding bands at a given (SSB) reference GPS-time .
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
Definition: SFTfileIO.h:238