Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
ppe_readdata.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2014 Matthew Pitkin
3*
4* This program is free software; you can redistribute it and/or modify
5* it under the terms of the GNU General Public License as published by
6* the Free Software Foundation; either version 2 of the License, or
7* (at your option) any later version.
8*
9* This program is distributed in the hope that it will be useful,
10* but WITHOUT ANY WARRANTY; without even the implied warranty of
11* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12* GNU General Public License for more details.
13*
14* You should have received a copy of the GNU General Public License
15* along with with program; see the file COPYING. If not, write to the
16* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17* MA 02110-1301 USA
18*/
19
20
21/******************************************************************************/
22/* DATA READING FUNCTIONS */
23/******************************************************************************/
24
25#include "config.h"
26#include "ppe_readdata.h"
27
28/**
29 * \brief Reads in the input gravitational wave data files, or creates fake data sets.
30 *
31 * The function will check whether data files are being input of fake data is to be generated. If using real data the \c
32 * detectors command line input must list the names of the detectors from which each of the data sets comes, with names
33 * separated by commas - allowed detector names are H1, H2, L1, V1, G1 or T1, for LIGO Hanford Observatory 1, LIGO
34 * Hanford Observatory 2, LIGO Livingston Observatory, Virgo, GEO600/HF and TAMA300 respectively. If requiring fake data
35 * to be generated then the \c fake-data command line argument must list the detectors for which fake data is required
36 * (again separated by commas) - these can be the same names as above, although if an 'A' is appended to the LIGO of
37 * Virgo detector name then the Advanced detector is assumed (for use if generating data from designed sensitivity
38 * curves).
39 *
40 * If generating fake data the noise floor can either be created using models of the detector design sensitivities (if
41 * just \c fake-data is set), or with noise levels defined using the \c fake-psd command line argument. If using \c
42 * fake-psd it should list the signle-sided power spectral density required for each detector listed in \c fake-data
43 * (again separated by commas). If using the design sensitivity models the \c par-file argument will be used to find the
44 * noise at the correct frequency, which is here assumed to be twice the rotation frequency. The start time (in GPS
45 * seconds), data length (in seconds) and sample interval (in seconds) can be specified for each fake data set by
46 * supplying comma separated lists to the \c fake-starts, \c fake-lengths and \c fake-dt command line arguments. By
47 * default these values are GPS 900000000 (13th July 2008 at 15:59:46 UTC), 86400 seconds (1 solar day) and 60 seconds
48 * (1/60 Hz sample rate) respectively. The fake data is drawn from a normal distribution using \c XLALNormalDeviates and
49 * scaled with the appropriate PSD.
50 *
51 * The number of data files required to be read in, or number of fake data sets generated depends on the pulsar model
52 * type, which is specified by the number of frequency harmonics given by the \c harmonics command line argument. This
53 * should be a list of comma separated values giving the frequency of the signal harmonics to be included. E.g.
54 * <ol>
55 * <li>For a model with \f$ l=m=2 \f$ (i.e. a triaxial star with a signal defined in e.g. \cite DupuisWoan2005), which
56 * purely emits at twice the rotation frequency, the \c harmonics input would just be \c 2. This requires that for each
57 * specified detector there is <b>one</b> input file containing data heterodyned at twice the rotation frequency of the
58 * pulsar.</li>
59 * <li>For a model including the two harmonics \f$ l=2 \f$ , \f$ m=1,2 \f$ (see e.g. \cite Jones2010), which produces
60 * emission at both the rotation frequency <i>and</i> twice the rotation frequency, the \c harmonics input would be \c
61 * 1,2. This requires that for each specified detector there are two input files containing data heterodyned and the
62 * rotation frequency <i>and</i> twice the rotation frequency (these must be in the same order as the harmonics
63 * are listed).</li>
64 * </ol>
65 * The default model, if none is specified, is the triaxial source model with emission at just twice the rotation
66 * frequency. At the moment only these two models above can be used, although this could be expanded in the future.
67 *
68 * If creating fake data for a specific model then the number of PSDs, start time, lengths and sample intervals
69 * specified must be equivalent to the number of input files that would have been required e.g. if using the pinned
70 * superfluid model and requiring data for H1 and L1 then four fake PSDs would be required (the first pair at the
71 * pulsars rotation frequency and twice that in H1, and the seconds pair at the pulsars rotation frequency and twice
72 * that in L1). These most be specified in the same order as the detectors.
73 *
74 * Any new models added can require and arbitrary number of inputs for a given detector, however the heterodyned
75 * frequencies of each input must be hardcoded into the function.
76 *
77 * If using real data the files must be specified in the \c input-files command line argument - these should be comma
78 * separated for multiple files and be in the same order at the associated detector from which they came given by the
79 * \c detectors command. Any potentially spuriously large data values, with an absolute value of greater than 1e-18
80 * will be vetoed when the data in read in. If requiring a different threshold for the veto (e.g., you are working
81 * with very low frequency data, or simulated data with a different scaling) then the value can be changed by setting
82 * it with the \c --veto-threshold command line argument.
83 *
84 * The function also checks that valid Earth and Sun ephemeris files (from the lalpulsar suite) are set with the \c
85 * ephem-earth and \c ephem-sun arguments, and that a valid output file for the nested samples is set via the \c outfile
86 * argument.
87 *
88 * The function will by default also call \c chop_n_merge() for each data set, which will split the data into chunks
89 * during which it can be considered Gaussian and stationary. The command line arguments \c chunk-min and \c chunk-max
90 * can be used to specify hardwire minimum and maximum lengths of chunks that are allowable. By default the maximum
91 * chunk length is 0, which corresponds to no maximum value being set. If the \c --oldChunks flag is set then data will
92 * be split as in the older version of the parameter estimation code, in which the chunk length is fixed, except for the
93 * possibility of shorter segments at the end of science segments.
94 *
95 * The \c log_factorial array is also filled in with values of the log of the factorial of all numbers up to the maximum
96 * length of the data. This is used in likelihood calculations.
97 *
98 * \param runState [in] A pointer to the LALInferenceRunState
99 */
101{
102 ProcessParamsTable *ppt = NULL, *ppt2 = NULL;
103 ProcessParamsTable *commandLine = runState->commandLine;
104
105 CHAR *detectors = NULL;
106 CHAR *inputfile = NULL;
107
108 CHAR *filestr = NULL;
109
110 CHAR *efile = NULL, *sfile = NULL, *tfile = NULL;
111 TimeCorrectionType ttype;
112
113 CHAR *tempdets = NULL;
114 CHAR *tempdet = NULL;
115
116 REAL8 *fpsds = NULL;
117 CHAR *fakestarts = NULL, *fakelengths = NULL, *fakedt = NULL;
118 REAL8 *fstarts = NULL, *flengths = NULL, *fdt = NULL;
119
120 CHAR dets[MAXDETS][256];
121 INT4 numDets = 1, i = 0, numPsds = 0, numLengths = 0, numStarts = 0;
122 INT4 numDt = 0, count = 0;
123 UINT4 maxlen = 0;
124
125 LALInferenceIFOData *ifodata = NULL;
126 LALInferenceIFOData *prev = NULL;
127
128 LALInferenceIFOModel *ifomodel = NULL;
129 LALInferenceIFOModel *prevmodel = NULL;
130
131 UINT4 seed = 0; /* seed for data generation */
132 RandomParams *randomParams = NULL;
133
134 CHAR *harmonics = NULL;
135 REAL8Vector *modelFreqFactors = NULL;
136 INT4 ml = 1;
137
138 CHAR *parFile = NULL;
140
141 runState->data = NULL;
142
143 /* Initialize the model, as it will hold IFO params and signal buffers */
144 /* single thread */
145 runState->threads[0].model = XLALCalloc( 1, sizeof( LALInferenceModel ) );
146
147 /* timing values */
148 struct timeval time1, time2;
149 REAL8 tottime;
150
151 if ( LALInferenceCheckVariable( runState->algorithmParams, "timefile" ) ) {
152 gettimeofday( &time1, NULL );
153 }
154
155 /* check pulsar model required by getting the frequency harmonics */
156 ppt = LALInferenceGetProcParamVal( commandLine, "--harmonics" );
157 if ( ppt ) {
159 } else {
160 harmonics = XLALStringDuplicate( "2" ); /* default to using harmonic at twice the rotation rate */
161 }
162
163 CHAR *tmpharms = NULL, *tmpharm = NULL, harmval[256];
164
165 ml = count_csv( harmonics );
166 modelFreqFactors = XLALCreateREAL8Vector( ml );
167
168 tmpharms = XLALStringDuplicate( harmonics );
169 for ( i = 0; i < ml; i++ ) {
170 tmpharm = XLALStringToken( &tmpharms, ",", 0 );
171 XLALStringCopy( harmval, tmpharm, strlen( tmpharm ) + 1 );
172
173 modelFreqFactors->data[i] = atof( harmval );
174 }
175 XLALFree( tmpharms );
177
178 ppt = LALInferenceGetProcParamVal( runState->commandLine, "--par-file" );
179 if ( ppt == NULL ) {
180 fprintf( stderr, "Must specify --par-file!\n" );
181 exit( 1 );
182 }
183 parFile = ppt->value;
184
185 /* get the pulsar parameters to give a value of f */
186 pulsar = XLALReadTEMPOParFile( parFile );
187
188 /* get the detectors - must */
189 ppt = LALInferenceGetProcParamVal( commandLine, "--detectors" );
190 ppt2 = LALInferenceGetProcParamVal( commandLine, "--fake-data" );
191 if ( ppt && !ppt2 ) {
192 detectors = ppt->value;
193
194 /* count the number of detectors from command line argument of comma separated vales and set their names */
195 tempdets = XLALStringDuplicate( detectors );
196
197 if ( ( numDets = count_csv( detectors ) ) > MAXDETS ) {
198 fprintf( stderr, "Error... too many detectors specified. Increase MAXDETS to be greater than %d if necessary.\n", MAXDETS );
199 exit( 0 );
200 }
201
202 for ( i = 0; i < numDets; i++ ) {
203 tempdet = XLALStringToken( &tempdets, ",", 0 );
204 XLALStringCopy( dets[i], tempdet, strlen( tempdet ) + 1 );
205 }
206 }
207 /*Get psd values for generating fake data.*/
208 /*=========================================================================*/
209 /*if using fake data and no detectors are specified*/
210 else if ( ppt2 && !ppt ) {
211 detectors = ppt2->value;
212
213 fpsds = XLALCalloc( MAXDETS * ml, sizeof( REAL8 ) );
214
215 if ( ( numDets = count_csv( detectors ) ) > MAXDETS ) {
216 fprintf( stderr, "Error... too many detectors specified. Increase MAXDETS to be greater than %d if necessary.\n", MAXDETS );
217 exit( 0 );
218 }
219
220 /*------------------------------------------------------------------------*/
221 /* get noise psds if specified */
222 ppt = LALInferenceGetProcParamVal( commandLine, "--fake-psd" );
223 if ( ppt ) {
224 CHAR *psds = NULL, *tmppsds = NULL, *tmppsd = NULL, psdval[256];
225
226 psds = ppt->value;
227
228 tmppsds = XLALStringDuplicate( psds );
229 tempdets = XLALStringDuplicate( detectors );
230
231 /* count the number of PSDs (comma seperated values) to compare to number of detectors */
232 if ( ( numPsds = count_csv( psds ) ) != ml * numDets ) {
233 fprintf( stderr, "Error... for %d harmonics the number of PSDs for fake data must be %d times the number of \
234detectors specified (no. dets = %d)\n", ml, ml, numDets );
235 exit( 0 );
236 }
237
238 for ( i = 0; i < ml * numDets; i++ ) {
239 CHAR *tmpstr = NULL;
240
241 tmppsd = XLALStringToken( &tmppsds, ",", 0 );
242 XLALStringCopy( psdval, tmppsd, strlen( tmppsd ) + 1 );
243 fpsds[i] = atof( psdval );
244
245 /* set detector */
246 if ( i % ml == 0 ) {
247 tempdet = XLALStringToken( &tempdets, ",", 0 );
248
249 if ( ( tmpstr = strstr( tempdet, "A" ) ) != NULL ) { /* have advanced */
250 XLALStringCopy( dets[FACTOR( i, ml )], tmpstr + 1, strlen( tmpstr ) + 1 );
251 } else {
252 XLALStringCopy( dets[FACTOR( i, ml )], tempdet, strlen( tempdet ) + 1 );
253 }
254 }
255 }
256 }
257 /*------------------------------------------------------------------------*/
258 else { /* get PSDs from model functions and set detectors */
259 REAL8 pfreq = 0.;
260
261 /* putting in pulsar frequency at f here */
262 if ( PulsarCheckParam( pulsar, "F" ) ) {
264 } else {
265 XLAL_ERROR_VOID( XLAL_EINVAL, "No source frequency given in parameter file" );
266 }
267
268 tempdets = XLALStringDuplicate( detectors );
269
270 for ( i = 0; i < ml * numDets; i++ ) {
271 CHAR *tmpstr = NULL;
272 REAL8 psdvalf = 0.;
273
274 numPsds++;
275
276 if ( i % ml == 0 ) {
277 tempdet = XLALStringToken( &tempdets, ",", 0 );
278 }
279
280 if ( ( tmpstr = strstr( tempdet, "A" ) ) != NULL ) { /* have Advanced */
281 XLALStringCopy( dets[FACTOR( i, ml )], tmpstr + 1, strlen( tmpstr ) + 1 );
282
283 if ( !strcmp( dets[FACTOR( i, ml )], "H1" ) || !strcmp( dets[FACTOR( i, ml )], "L1" ) || !strcmp( dets[FACTOR( i, ml )], "H2" ) ) { /* ALIGO */
284 psdvalf = XLALSimNoisePSDaLIGOZeroDetHighPower( pfreq * modelFreqFactors->data[( INT4 )( fmod( i, ml ) )] );
285 } else if ( !strcmp( dets[FACTOR( i, ml )], "V1" ) ) { /* AVirgo */
286 psdvalf = XLALSimNoisePSDAdvVirgo( pfreq * modelFreqFactors->data[( INT4 )( fmod( i, ml ) )] );
287 } else {
288 fprintf( stderr, "Error... trying to use Advanced detector that is not available!\n" );
289 exit( 0 );
290 }
291 } else { /* initial detector */
292 XLALStringCopy( dets[FACTOR( i, ml )], tempdet, strlen( tempdet ) + 1 );
293
294 if ( !strcmp( dets[FACTOR( i, ml )], "H1" ) || !strcmp( dets[FACTOR( i, ml )], "L1" ) || !strcmp( dets[FACTOR( i, ml )], "H2" ) ) { /* Initial LIGO */
295 psdvalf = XLALSimNoisePSDiLIGOSRD( pfreq * modelFreqFactors->data[( INT4 )( fmod( i, ml ) )] );
296
297 /* divide H2 psds by 2 */
298 if ( !strcmp( dets[FACTOR( i, ml )], "H2" ) ) {
299 psdvalf /= 2.;
300 }
301 } else if ( !strcmp( dets[FACTOR( i, ml )], "V1" ) ) { /* Initial Virgo */
302 psdvalf = XLALSimNoisePSDVirgo( pfreq * modelFreqFactors->data[( INT4 )( fmod( i, ml ) )] );
303 } else if ( !strcmp( dets[FACTOR( i, ml )], "G1" ) ) { /* GEO 600 */
304 psdvalf = XLALSimNoisePSDGEO( pfreq * modelFreqFactors->data[( INT4 )( fmod( i, ml ) )] );
305 } else if ( !strcmp( dets[FACTOR( i, ml )], "T1" ) ) { /* TAMA300 */
306 psdvalf = XLALSimNoisePSDTAMA( pfreq * modelFreqFactors->data[( INT4 )( fmod( i, ml ) )] );
307 } else {
308 fprintf( stderr, "Error... trying to use detector that is not available!\n" );
309 exit( 0 );
310 }
311 }
312
313 fpsds[i] = psdvalf;
314 }
315 }
316 /*generate the fake data timestamps.*/
317 /*====================================================================*/
318
319 fstarts = XLALCalloc( MAXDETS * ml, sizeof( REAL8 ) );
320 ppt = LALInferenceGetProcParamVal( commandLine, "--fake-starts" );
321 if ( ppt ) {
322 CHAR *tmpstarts = NULL, *tmpstart = NULL, startval[256];
323 fakestarts = ppt->value;
324
325 if ( ( numStarts = count_csv( fakestarts ) ) != numDets * ml ) {
326 fprintf( stderr, "Error... for %d harmonics the number of start times for fake data must be %d times the number \
327of detectors specified (no. dets = %d)\n", ml, ml, numDets );
328 exit( 0 );
329 }
330
331 tmpstarts = XLALStringDuplicate( fakestarts );
332 for ( i = 0; i < ml * numDets; i++ ) {
333 tmpstart = XLALStringToken( &tmpstarts, ",", 0 );
334 XLALStringCopy( startval, tmpstart, strlen( tmpstart ) + 1 );
335
336 fstarts[i] = atof( startval );
337 }
338 } else { /* set default GPS 900000000 - 13th July 2008 at 15:59:46 */
339 for ( i = 0; i < ml * numDets; i++ ) {
340 fstarts[i] = 900000000.;
341 }
342 }
343
344 flengths = XLALCalloc( MAXDETS * ml, sizeof( REAL8 ) );
345 ppt = LALInferenceGetProcParamVal( commandLine, "--fake-lengths" );
346 if ( ppt ) {
347 CHAR *tmplengths = NULL, *tmplength = NULL, lengthval[256];
348 fakelengths = ppt->value;
349
350 if ( ( numLengths = count_csv( fakelengths ) ) != numDets * ml ) {
351 fprintf( stderr, "Error... for %d harmonics the number of data lengths for fake data must be %d times the \
352number of detectors specified (no. dets = %d)\n", ml, ml, numDets );
353 exit( 0 );
354 }
355
356 tmplengths = XLALStringDuplicate( fakelengths );
357 for ( i = 0; i < ml * numDets; i++ ) {
358 tmplength = XLALStringToken( &tmplengths, ",", 0 );
359 XLALStringCopy( lengthval, tmplength, strlen( tmplength ) + 1 );
360 flengths[i] = atof( lengthval );
361 }
362 } else { /* set default (86400 seconds or 1 day) */
363 for ( i = 0; i < ml * numDets; i++ ) {
364 flengths[i] = 86400.;
365 }
366 }
367
368 fdt = XLALCalloc( MAXDETS * ml, sizeof( REAL8 ) );
369 ppt = LALInferenceGetProcParamVal( commandLine, "--fake-dt" );
370 if ( ppt ) {
371 CHAR *tmpdts = NULL, *tmpdt = NULL, dtval[256];
372 fakedt = ppt->value;
373
374 if ( ( numDt = count_csv( fakedt ) ) != ml * numDets ) {
375 fprintf( stderr, "Error... for %d harmonics the number of sample time steps for fake data must be %d times the \
376number of detectors specified (no. dets =%d)\n", ml, ml, numDets );
377 exit( 0 );
378 }
379
380 tmpdts = XLALStringDuplicate( fakedt );
381
382 for ( i = 0; i < ml * numDets; i++ ) {
383 tmpdt = XLALStringToken( &tmpdts, ",", 0 );
384 XLALStringCopy( dtval, tmpdt, strlen( tmpdt ) + 1 );
385 fdt[i] = atof( dtval );
386 }
387 } else { /* set default (60 sesonds) */
388 for ( i = 0; i < ml * numDets; i++ ) {
389 fdt[i] = 60.;
390 }
391 }
392 }
393 /*psds set and timestamps set.*/
394 /*====================================================================*/
395 else {
396 fprintf( stderr, "Error... --detectors OR --fake-data needs to be set.\n" );
397 fprintf( stderr, USAGE, commandLine->program );
398 exit( 0 );
399 }
400
401 XLALFree( tempdets );
402
403 runState->threads[0].model->ifo_loglikelihoods = XLALMalloc( sizeof( REAL8 ) * ml * numDets );
404 runState->threads[0].model->ifo_SNRs = XLALMalloc( sizeof( REAL8 ) * ml * numDets );
405
406 UINT4 nstreams = ml * numDets;
408
409 /* check the output file is given */
410 if ( !LALInferenceGetProcParamVal( commandLine, "--outfile" ) ) {
411 fprintf( stderr, "Error... --outfile needs to be set.\n" );
412 fprintf( stderr, USAGE, commandLine->program );
413 exit( 0 );
414 }
415
416 ppt = LALInferenceGetProcParamVal( commandLine, "--input-files" );
417 if ( ppt ) {
418 inputfile = ppt->value;
419 }
420
421 if ( ( inputfile == NULL || strlen( inputfile ) == 0 ) && !ppt2 ) {
422 fprintf( stderr, "Error... an input file or fake data needs to be set.\n" );
423 fprintf( stderr, USAGE, commandLine->program );
424 exit( 0 );
425 }
426
427 /* count the number of input files (by counting commas) and check it's equal to twice the number of detectors */
428 if ( !ppt2 ) { /* if using real data */
429 count = count_csv( inputfile );
430
431 if ( count != ml * numDets ) {
432 fprintf( stderr, "Error... for %d harmonics the number of input files given must be %d times the number of \
433detectors specified (no. dets =%d)\n", ml, ml, numDets );
434 exit( 0 );
435 }
436 }
437
438 /* reset filestr if using real data (i.e. not fake) */
439 if ( !ppt2 ) {
440 filestr = XLALStringDuplicate( inputfile );
441 }
442
443 /* set random number generator in case when that fake data is used */
444 ppt = LALInferenceGetProcParamVal( commandLine, "--randomseed" );
445 if ( ppt != NULL ) {
446 seed = atoi( ppt->value );
447 } else {
448 seed = 0; /* will be set from system clock */
449 }
450
451 /* set flags for whether noise sigma is input or required to be calculated */
452 INT4 inputsigma = 0, computesigma = 0, gaussianLike = 0;
453
454 if ( LALInferenceGetProcParamVal( commandLine, "--gaussian-like" ) ) {
455 gaussianLike = 1;
456 }
457
458 /* set reheterodyne frequency (0 if no reheterodyne is required) */
459 REAL8 rehetfreq = 0;
460 if ( LALInferenceGetProcParamVal( commandLine, "--reheterodyne" ) ) {
461 if ( *LALInferenceGetProcParamVal( commandLine, "--reheterodyne" )->value == '\0' ) {
462 fprintf( stderr, "Error... --reheterodyne needs frequency as argument.\n" );
463 fprintf( stderr, "Provide argument or remove flag\n." );
464 exit( 1 );
465 } else {
466 rehetfreq = atof( LALInferenceGetProcParamVal( commandLine, "--reheterodyne" )->value );
467 }
468 }
469
470 /* ignore data with absolute value larger than some input value (defaults to 1e-18) to chop out spurious extremely load noise */
471 ppt = LALInferenceGetProcParamVal( commandLine, "--veto-threshold" );
472 REAL8 vetothresh = 1e-18;
473 if ( ppt != NULL ) {
474 vetothresh = atof( ppt->value );
475 }
476
477 /* initialise random number generator if seed is zero */
478 if ( seed == 0 ) {
479 randomParams = XLALCreateRandomParams( 0 );
480 }
481
482 /* check if truncating data or only using part of it - the "--start-time" and
483 "--end-time" flags can be used to specify a chunk of data to use, and these
484 will overrule any "--truncate-*" flags that are used. This flags will only
485 be used if using real data, but not if generating fake data. */
486 REAL8 startTimeValue = 0., endTimeValue = INFINITY;
487 if ( LALInferenceGetProcParamVal( commandLine, "--start-time" ) ) {
488 startTimeValue = atof( LALInferenceGetProcParamVal( commandLine, "--start-time" )->value );
489 }
490 if ( LALInferenceGetProcParamVal( commandLine, "--end-time" ) ) {
491 endTimeValue = atof( LALInferenceGetProcParamVal( commandLine, "--end-time" )->value );
492 }
493
494 if ( endTimeValue <= startTimeValue ) {
495 fprintf( stderr, "Error... start time is before end time.\n" );
496 exit( 1 );
497 }
498
499 if ( ( startTimeValue == 0 ) && ( endTimeValue == INFINITY ) ) {
500 /* there are three truncation modes: "--truncate-time" uses a maximum */
501 /* GPS time, "--truncate-samples" uses maximum number of samples and */
502 /* "--truncate-fraction" uses a fraction of the number of samples. */
503 UINT4 truncate = 0;
504 char truncationFlags[][25] = { "--truncate-time", "--truncate-samples", "--truncate-fraction" };
505
506 size_t trunci = 0;
507 for ( trunci = 0; trunci < sizeof( truncationFlags ) / sizeof( truncationFlags[0] ); trunci++ ) {
508 if ( LALInferenceGetProcParamVal( commandLine, truncationFlags[trunci] ) ) {
509 if ( *LALInferenceGetProcParamVal( commandLine, truncationFlags[trunci] )-> value != '\0' ) {
510 truncate++;
511 endTimeValue = atof( LALInferenceGetProcParamVal( commandLine, truncationFlags[trunci] )->value );
512 } else {
513 fprintf( stderr, "Error... truncation option requires a value.\n" );
514 fprintf( stderr, USAGE, commandLine->program );
515 exit( 1 );
516 }
517 }
518 }
519
520 if ( truncate > 1 ) {
521 fprintf( stderr, "Error... can only take one truncation flag.\n" );
522 fprintf( stderr, USAGE, commandLine->program );
523 exit( 1 );
524 }
525 }
526
527 /* read in data, needs to read in two sets of data for each ifo for pinsf model */
528 for ( i = 0, prev = NULL, prevmodel = NULL ; i < ml * numDets ; i++, prev = ifodata, prevmodel = ifomodel ) {
529 CHAR *datafile = NULL;
530 REAL8 times = 0;
532 REAL8 dataValsRe = 0., dataValsIm = 0., sigmaVals = 0.;
533 REAL8Vector *temptimes = NULL;
534 UINT4 j = 0, k = 0, datalength = 0;
535 ProcessParamsTable *ppte = NULL, *ppts = NULL, *pptt = NULL;
536
537 CHAR *filebuf = NULL;
538
539 count = 0;
540
541 /* initialise random number generator (if using a non-zero seed).
542 * Moved into det loop so same random seed can be used with
543 * different detector combos and still get same noise realisation */
544 if ( seed != 0 ) {
545 randomParams = XLALCreateRandomParams( seed + i );
546 }
547
548 ifodata = XLALCalloc( 1, sizeof( LALInferenceIFOData ) );
549 ifodata->likeli_counter = 0;
550 ifodata->templa_counter = 0;
551 ifodata->next = NULL;
552
553 ifomodel = XLALMalloc( sizeof( LALInferenceIFOModel ) );
554 ifomodel->extraData = XLALCalloc( 1, sizeof( IFOModelExtraData ) );
555 ifomodel->params = XLALCalloc( 1, sizeof( LALInferenceVariables ) );
556 ifomodel->next = NULL;
557
558 /* add frequency factors variable */
559 REAL8Vector *freqFactorsCopy = NULL;
560 freqFactorsCopy = XLALCreateREAL8Vector( modelFreqFactors->length );
561 memcpy( freqFactorsCopy->data, modelFreqFactors->data, sizeof( REAL8 )*modelFreqFactors->length );
562 LALInferenceAddVariable( ifomodel->params, "freqfactors", &freqFactorsCopy, LALINFERENCE_REAL8Vector_t, LALINFERENCE_PARAM_FIXED );
563
564 /* check if using non-GR model */
565 if ( LALInferenceGetProcParamVal( commandLine, "--nonGR" ) ) {
566 UINT4 nonGR = 1;
568 /* check if nonGR argument has an associated input variable */
569 if ( *LALInferenceGetProcParamVal( commandLine, "--nonGR" )->value != '\0' ) {
570 CHAR *nonGRmodel = NULL;
571 nonGRmodel = XLALStringDuplicate( LALInferenceGetProcParamVal( commandLine, "--nonGR" )->value );
572 LALInferenceAddVariable( ifomodel->params, "nonGRmodel", &nonGRmodel, LALINFERENCE_string_t, LALINFERENCE_PARAM_FIXED );
573 }
574 }
575
576 if ( i == 0 ) {
577 runState->data = ifodata;
578 runState->threads[0].model->ifo = ifomodel;
579 }
580 if ( i > 0 ) {
581 prev->next = ifodata;
582 prevmodel->next = ifomodel;
583 }
584
585 /* set detector */
586 ifodata->detector = XLALCalloc( 1, sizeof( *ifodata->detector ) );
587 memcpy( ifodata->detector, XLALGetSiteInfo( dets[FACTOR( i, ml )] ), sizeof( *ifodata->detector ) );
588 ifomodel->detector = XLALCalloc( 1, sizeof( *ifomodel->detector ) );
589 memcpy( ifomodel->detector, XLALGetSiteInfo( dets[FACTOR( i, ml )] ), sizeof( *ifomodel->detector ) );
590 strncpy( ifodata->name, dets[FACTOR( i, ml )], DETNAMELEN - 1 );
591
592 /* set dummy initial time */
593 gpstime.gpsSeconds = 0;
594 gpstime.gpsNanoSeconds = 0;
595
596 /* allocate time domain data - will be dynamically allocated as data read*/
597 ifodata->compTimeData = NULL;
598 ifodata->compTimeData = XLALCreateCOMPLEX16TimeSeries( "", &gpstime, 0., 1., &lalSecondUnit, 1 );
599
600 /* always allocate memory for variances (as these will always be needed in calculating SNRs) */
601 ifodata->varTimeData = NULL;
602 ifodata->varTimeData = XLALCreateREAL8TimeSeries( "", &gpstime, 0., 1., &lalSecondUnit, 1 );
603
604 /* allocate time domain model */
605 ifomodel->compTimeSignal = NULL;
606 ifomodel->compTimeSignal = XLALCreateCOMPLEX16TimeSeries( "", &gpstime, 0., 1., &lalSecondUnit, 1 );
607
608 /*============================ GET DATA ==================================*/
609 /* get i'th filename from the comma separated list */
610 if ( !ppt2 ) { /* if using real data read in from the file */
611 datafile = XLALStringToken( &filestr, ",", 0 );
612
613 j = 0;
614
615 /* read in data */
616 temptimes = XLALCreateREAL8Vector( 1 );
617
618 /* read in all the data then ignore lines starting with # or % */
619 filebuf = XLALFileLoad( datafile );
620
621 /* separate data into lines */
622 TokenList *tlist = NULL;
623 if ( XLALCreateTokenList( &tlist, filebuf, "\n" ) != XLAL_SUCCESS ) {
624 fprintf( stderr, "Error... could not convert data into separate lines.\n" );
625 exit( 3 );
626 }
627
628 INT4 nvals = 0; /* number of values in a line */
629
630 for ( k = 0; k < tlist->nTokens; k++ ) {
631 /* search for a comment character in the string */
632 if ( strchr( tlist->tokens[k], '#' ) || strchr( tlist->tokens[k], '%' ) ) {
633 continue;
634 } else { /* read in data from string */
635 if ( j == 0 ) {
636 /* check the number of values in the line by counting the number of value separated by whitespace */
637 TokenList *tline = NULL;
638 XLALCreateTokenList( &tline, tlist->tokens[k], " \t" );
639 nvals = ( INT4 )tline->nTokens;
640 XLALDestroyTokenList( tline );
641
642 /* set whether we need to compute variance, or whether the standard deviation is in the input file */
643 if ( nvals == 3 && gaussianLike ) {
644 computesigma = 1;
645 }
646 if ( nvals == 4 ) {
647 inputsigma = 1;
648 }
649 }
650
651 if ( nvals == 3 ) { /* no sigma value is in the input file */
652 int rc = sscanf( tlist->tokens[k], "%lf%lf%lf", &times, &dataValsRe, &dataValsIm );
653 if ( rc != nvals ) {
654 continue; /* ignore the line */
655 }
656 } else if ( nvals == 4 ) { /* sigma is in the input file */
657 int rc = sscanf( tlist->tokens[k], "%lf%lf%lf%lf", &times, &dataValsRe, &dataValsIm, &sigmaVals );
658 if ( rc != nvals ) {
659 continue; /* ignore the line */
660 }
661 } else {
662 fprintf( stderr, "Error... unrecognised number of values in first line of data file %s.\n", datafile );
663 exit( 3 );
664 }
665 /* ignore excessively large spurious values as they can screw things up */
666 if ( fabs( dataValsRe ) > vetothresh || fabs( dataValsIm ) > vetothresh ) {
667 continue;
668 }
669
670 /* make sure timestamps are unique */
671 if ( j > 0 ) {
672 UINT4 notunique = 0;
673 for ( UINT4 ucount = 0; ucount < j; ucount++ ) {
674 if ( times == temptimes->data[ucount] ) {
675 notunique = 1;
676 break;
677 }
678 }
679 if ( notunique ) {
680 continue;
681 }
682 }
683 }
684
685 /* check whether to only inlcude a section of the data */
686 if ( !isinf( endTimeValue ) || ( startTimeValue != 0. ) ) {
687 if ( times < startTimeValue ) {
688 continue;
689 }
690
691 if ( LALInferenceGetProcParamVal( commandLine, "--end-time" ) ||
692 LALInferenceGetProcParamVal( commandLine, "--truncate-time" ) ) {
693 /* assume endTimeValue is a GPS time */
694 if ( times > endTimeValue ) {
695 /* exit loop */
696 break;
697 }
698 }
699
700 if ( LALInferenceGetProcParamVal( commandLine, "--truncate-samples" ) ) {
701 /* assume endTimeValue is an index */
702 if ( j >= ( UINT4 )endTimeValue ) {
703 /* exit loop */
704 break;
705 }
706 }
707 }
708
709 j++;
710
711 /* dynamically allocate more memory */
712 ifodata->compTimeData = XLALResizeCOMPLEX16TimeSeries( ifodata->compTimeData, 0, j );
714 ifodata->varTimeData = XLALResizeREAL8TimeSeries( ifodata->varTimeData, 0, j );
715
716 temptimes = XLALResizeREAL8Vector( temptimes, j );
717
718 /* Note: j-1 because we already added to j above */
719 temptimes->data[j - 1] = times;
720
721 /* reheterodyne data if required */
722 if ( rehetfreq != 0 ) {
723 /* create template */
724 REAL8 deltaPhase = 2. * LAL_PI * rehetfreq * times;
725 COMPLEX16 dataTemp = ( dataValsRe * cos( -deltaPhase ) - dataValsIm * sin( -deltaPhase ) ) +
726 I * ( dataValsRe * sin( -deltaPhase ) + dataValsIm * cos( -deltaPhase ) );
727 ifodata->compTimeData->data->data[j - 1] = dataTemp;
728 } else {
729 ifodata->compTimeData->data->data[j - 1] = dataValsRe + I * dataValsIm;
730 }
731
732 if ( inputsigma ) {
733 ifodata->varTimeData->data->data[j - 1] = SQUARE( sigmaVals );
734 }
735 }
736
737 if ( j == 0 ) {
738 fprintf( stderr, "Error... nothing read in from data file %s.\n", datafile );
739 exit( 3 );
740 }
741
742 XLALDestroyTokenList( tlist );
743 XLALFree( filebuf );
744
745 datalength = j;
746
747 /* truncate data if --truncate-fraction was passed */
748 if ( LALInferenceGetProcParamVal( commandLine, "--truncate-fraction" ) &&
749 !LALInferenceGetProcParamVal( commandLine, "--start-time" ) &&
750 !LALInferenceGetProcParamVal( commandLine, "--end-time" ) ) {
751 /* assume endTimeValue is a decimal between 0 and 1 */
752 /* (fraction of the number of samples) */
753 if ( endTimeValue <= 0 || endTimeValue > 1 ) {
754 fprintf( stderr, "Error... truncation fraction must be between 0 and 1" );
755 exit( 3 );
756 } else {
757 /* discard excess data */
758 int truncationIndex = floor( endTimeValue * datalength );
759 ifodata->compTimeData = XLALResizeCOMPLEX16TimeSeries( ifodata->compTimeData, 0, truncationIndex );
760 ifomodel->compTimeSignal = XLALResizeCOMPLEX16TimeSeries( ifomodel->compTimeSignal, 0, truncationIndex );
761 ifodata->varTimeData = XLALResizeREAL8TimeSeries( ifodata->varTimeData, 0, truncationIndex );
762 }
763 }
764
765 datalength = ifodata->compTimeData->data->length;
766
767 /* allocate data time stamps */
768 IFO_XTRA_DATA( ifomodel )->times = NULL;
769 IFO_XTRA_DATA( ifomodel )->times = XLALCreateTimestampVector( datalength );
770
771 /* fill in time stamps as LIGO Time GPS Vector */
772 for ( k = 0; k < datalength; k++ ) {
773 XLALGPSSetREAL8( &IFO_XTRA_DATA( ifomodel )->times->data[k], temptimes->data[k] );
774 }
775
776 /* sort temporary time vector into ascending order (to get minimum time difference) */
777 gsl_sort( temptimes->data, 1, datalength );
778
779 REAL8 sampledt = temptimes->data[1] - temptimes->data[0]; /* sample interval */
780 for ( k = 2; k < datalength; k++ ) {
781 if ( temptimes->data[i] - temptimes->data[i - 1] < sampledt ) {
782 sampledt = temptimes->data[i] - temptimes->data[i - 1];
783 }
784 }
785
786 XLALGPSSetREAL8( &ifodata->compTimeData->epoch, temptimes->data[0] );
787 XLALGPSSetREAL8( &ifomodel->compTimeSignal->epoch, temptimes->data[0] );
788 XLALGPSSetREAL8( &ifodata->varTimeData->epoch, temptimes->data[0] );
789
790 /* check whether to randomise the data by shuffling the time stamps (this will preserve the order of
791 * the data for working out stationary chunk, but randomise the signal) */
792 if ( LALInferenceGetProcParamVal( commandLine, "--randomise" ) ) {
793 INT4 randshufseed = atoi( LALInferenceGetProcParamVal( commandLine, "--randomise" )->value );
794 INT4 prevseed = gsl_rng_get( runState->GSLrandom ); // get previous RNG value
795 gsl_rng_set( runState->GSLrandom, randshufseed ); // set to value from randomise
796 gsl_ran_shuffle( runState->GSLrandom, &IFO_XTRA_DATA( ifomodel )->times->data[0], ( size_t )datalength, sizeof( LIGOTimeGPS ) ); // shuffle data times
797 gsl_rng_set( runState->GSLrandom, prevseed ); // reset to previous value
798 }
799
800 /* add data sample interval */
802
803 XLALDestroyREAL8Vector( temptimes );
804 } else { /* set up fake data */
805 /* if a Gaussian likelihood is required compute sigma from the data (to mimic real life) */
806 if ( gaussianLike ) {
807 computesigma = 1;
808 }
809
810 datalength = flengths[i] / fdt[i];
811
812 /* temporary real and imaginary data vectors */
813 REAL4Vector *realdata = NULL;
814 REAL4Vector *imagdata = NULL;
815
816 REAL8 psdscale = 0.;
817
818 /* allocate data time stamps */
819 IFO_XTRA_DATA( ifomodel )->times = NULL;
820 IFO_XTRA_DATA( ifomodel )->times = XLALCreateTimestampVector( ( UINT4 )datalength );
821
822 /* add data sample interval */
824
825 /* resize the data and model times series */
826 ifodata->compTimeData = XLALResizeCOMPLEX16TimeSeries( ifodata->compTimeData, 0, datalength );
827 ifomodel->compTimeSignal = XLALResizeCOMPLEX16TimeSeries( ifomodel->compTimeSignal, 0, datalength );
828 ifodata->varTimeData = XLALResizeREAL8TimeSeries( ifodata->varTimeData, 0, datalength );
829
830 /* create data drawn from normal distribution with zero mean and unit variance */
831 realdata = XLALCreateREAL4Vector( datalength );
832 imagdata = XLALCreateREAL4Vector( datalength );
833
834 XLALNormalDeviates( realdata, randomParams );
835 XLALNormalDeviates( imagdata, randomParams );
836
837 /* converts single sided psd into double sided psd, and then into a time domain noise standard deviation */
838 psdscale = sqrt( ( fpsds[i] / 2. ) / ( 2. * fdt[i] ) );
839
840 /* create time stamps and scale data with the PSD */
841 for ( k = 0; k < datalength; k++ ) {
842 /* set time stamp */
843 XLALGPSSetREAL8( &IFO_XTRA_DATA( ifomodel )->times->data[k], fstarts[i] + fdt[i] * ( REAL8 )k );
844
845 ifodata->compTimeData->data->data[k] = ( REAL8 )realdata->data[k] * psdscale + I * ( REAL8 )imagdata->data[k] * psdscale;
846 }
847
848 ifodata->compTimeData->epoch = IFO_XTRA_DATA( ifomodel )->times->data[0];
849 ifomodel->compTimeSignal->epoch = IFO_XTRA_DATA( ifomodel )->times->data[0];
850 ifodata->varTimeData->epoch = IFO_XTRA_DATA( ifomodel )->times->data[0];
851
852 XLALDestroyREAL4Vector( realdata );
853 XLALDestroyREAL4Vector( imagdata );
854 }
855
856 /* set ephemeris data */
857 IFO_XTRA_DATA( ifomodel )->ephem = XLALMalloc( sizeof( EphemerisData ) );
858
859 /* get ephemeris files */
860 ppte = LALInferenceGetProcParamVal( commandLine, "--ephem-earth" );
861 ppts = LALInferenceGetProcParamVal( commandLine, "--ephem-sun" );
862 pptt = LALInferenceGetProcParamVal( commandLine, "--ephem-timecorr" );
863
864 if ( ppte && ppts ) {
865 efile = XLALStringDuplicate( ppte->value );
866 sfile = XLALStringDuplicate( ppts->value );
867
868 if ( pptt ) {
869 tfile = XLALStringDuplicate( pptt->value );
870
871 if ( PulsarCheckParam( pulsar, "UNITS" ) ) {
872 if ( !strcmp( PulsarGetStringParam( pulsar, "UNITS" ), "TDB" ) ) {
873 ttype = TIMECORRECTION_TDB;
874 } else {
875 ttype = TIMECORRECTION_TCB; /* default to TCB otherwise */
876 }
877 } else {
878 ttype = TIMECORRECTION_TCB;
879 }
880 } else {
881 tfile = NULL;
883 }
884
885 /* check ephemeris files exist and if not output an error message */
886 if ( fopen( sfile, "r" ) == NULL || fopen( efile, "r" ) == NULL ) {
887 fprintf( stderr, "Error... ephemeris files not, or incorrectly, defined!\n" );
888 exit( 3 );
889 }
890 } else { /* try getting files automatically */
891 if ( !( ttype = XLALAutoSetEphemerisFiles( &efile, &sfile, &tfile, pulsar,
892 IFO_XTRA_DATA( ifomodel )->times->data[0].gpsSeconds, IFO_XTRA_DATA( ifomodel )->times->data[datalength - 1].gpsSeconds ) ) ) {
893 fprintf( stderr, "Error... not been able to set ephemeris files!\n" );
894 exit( 3 );
895 }
896 }
897
898 /* set up ephemeris information */
899 XLAL_CHECK_VOID( ( IFO_XTRA_DATA( ifomodel )->ephem = XLALInitBarycenter( efile, sfile ) ) != NULL, XLAL_EFUNC );
900 if ( tfile ) {
902 } else {
903 IFO_XTRA_DATA( ifomodel )->tdat = NULL;
904 }
905 IFO_XTRA_DATA( ifomodel )->ttype = ttype;
906
907 if ( efile ) {
908 XLALFree( efile );
909 }
910 if ( sfile ) {
911 XLALFree( sfile );
912 }
913 if ( tfile ) {
914 XLALFree( tfile );
915 }
916
917 if ( seed != 0 ) {
918 XLALDestroyRandomParams( randomParams );
919 }
920
921 /* get maximum data length */
922 if ( ifodata->compTimeData->data->length > maxlen ) {
923 maxlen = ifodata->compTimeData->data->length;
924 }
925 }
926
927 if ( seed == 0 ) {
928 XLALDestroyRandomParams( randomParams );
929 }
930 XLALFree( filestr );
931
932 /* chop the data into stationary chunks and also calculate the noise variance if required
933 * (note that if there is going to be a signal injected then this variance will be recalculated
934 * after the injection has been made to make the analysis most similar to a real case). */
935 UINT4 chunkMin, chunkMax;
936
937 /* Get chunk min and chunk max */
938 ppt = LALInferenceGetProcParamVal( commandLine, "--chunk-min" );
939 if ( ppt ) {
940 chunkMin = atoi( ppt->value );
941 } else {
942 chunkMin = CHUNKMIN; /* default minimum chunk length */
943 }
944
945 ppt = LALInferenceGetProcParamVal( commandLine, "--chunk-max" );
946 if ( ppt ) {
947 chunkMax = atoi( ppt->value );
948 } else {
949 chunkMax = CHUNKMAX; /* default maximum chunk length */
950 }
951
952 LALInferenceIFOData *datatmp = runState->data;
953 LALInferenceIFOModel *modeltmp = runState->threads[0].model->ifo;
954 while ( modeltmp ) {
955 UINT4Vector *chunkLength = NULL;
956
959
960 ppt = LALInferenceGetProcParamVal( commandLine, "--oldChunks" );
961 if ( ppt ) { /* use old style quasi-fixed data chunk lengths */
962 /* if a chunk max wasn't set use 30 mins by default */
963 if ( !LALInferenceGetProcParamVal( commandLine, "--chunk-max" ) ) {
964 chunkMax = 30;
965 LALInferenceSetVariable( modeltmp->params, "chunkMax", &chunkMax );
966 }
967
968 /* if sigma's have been input then there is just one chunk with a length of the full dataset */
969 if ( inputsigma ) {
970 chunkLength = XLALCreateUINT4Vector( 1 );
971 chunkLength->data[0] = datatmp->varTimeData->data->length;
972 } else {
973 chunkLength = get_chunk_lengths( modeltmp, chunkMax );
974 }
975 }
976 /* use new change points analysis to get chunks */
977 else {
978 /* if sigma's have been input then there is just one chunk with a length of the full dataset */
979 if ( inputsigma ) {
980 chunkLength = XLALCreateUINT4Vector( 1 );
981 chunkLength->data[0] = datatmp->varTimeData->data->length;
982 } else {
983 UINT4 outputchunks = 0;
984 if ( LALInferenceGetProcParamVal( commandLine, "--output-chunks" ) ) {
985 outputchunks = 1;
986 }
987 chunkLength = chop_n_merge( datatmp, chunkMin, chunkMax, outputchunks );
988 }
989 }
990
991 LALInferenceAddVariable( modeltmp->params, "chunkLength", &chunkLength, LALINFERENCE_UINT4Vector_t, LALINFERENCE_PARAM_FIXED );
992 LALInferenceAddVariable( modeltmp->params, "inputSigma", &inputsigma, LALINFERENCE_INT4_t, LALINFERENCE_PARAM_FIXED );
993
994 /* compute the variance */
995 if ( !inputsigma ) {
996 compute_variance( datatmp, modeltmp );
997 }
998
999 /* if we have a Gaussian likelihood and have computed sigma now we can reset the
1000 * number of chunks to be just one chunk (so we only need to sum over one data chunk) */
1001 if ( computesigma ) {
1002 UINT4Vector *chunkLengthNew = NULL;
1003 LALInferenceRemoveVariable( modeltmp->params, "chunkLength" );
1004
1005 chunkLengthNew = XLALCreateUINT4Vector( 1 );
1006 chunkLengthNew->data[0] = datatmp->varTimeData->data->length;
1007
1008 LALInferenceAddVariable( modeltmp->params, "chunkLength", &chunkLengthNew, LALINFERENCE_UINT4Vector_t, LALINFERENCE_PARAM_FIXED );
1009 }
1010
1011 /* set whether using Gaussian likelihood */
1012 if ( gaussianLike ) {
1013 LALInferenceAddVariable( modeltmp->params, "gaussianLikelihood", &gaussianLike, LALINFERENCE_INT4_t, LALINFERENCE_PARAM_FIXED );
1014 }
1015
1016 datatmp = datatmp->next;
1017 modeltmp = modeltmp->next;
1018 }
1019
1020 /* free memory */
1021 XLALFree( fdt );
1022 XLALFree( flengths );
1023 XLALFree( fstarts );
1024 XLALFree( fpsds );
1025
1027
1028 if ( LALInferenceCheckVariable( runState->algorithmParams, "timefile" ) ) {
1029 gettimeofday( &time2, NULL );
1030
1031 FILE *timefile = *( FILE ** )LALInferenceGetVariable( runState->algorithmParams, "timefile" );
1032 UINT4 timenum = *( UINT4 * )LALInferenceGetVariable( runState->algorithmParams, "timenum" );
1033 tottime = ( REAL8 )( ( time2.tv_sec + time2.tv_usec * 1.e-6 ) - ( time1.tv_sec + time1.tv_usec * 1.e-6 ) );
1034 fprintf( timefile, "[%d] %s: %.9le secs\n", timenum, __func__, tottime );
1035 timenum++;
1036 check_and_add_fixed_variable( runState->algorithmParams, "timenum", &timenum, LALINFERENCE_UINT4_t );
1037 }
1038}
1039
1040
1041/**
1042 * \brief Reads in the parameters of the pulsar being searched for
1043 *
1044 * This function reads in a pulsars parameters from the specified TEMPO-style .par file given by \c par-file using \c
1045 * XLALReadTEMPOParFile. This file must be specified and should contain at least the pulsars frequency, right
1046 * ascension and declination (any value not included will be zero by default). The file should contain the parameters
1047 * with which the detector data was heterodyned, as these are used to produce a signal phase template based on this
1048 * assumption.
1049 *
1050 * A example .par file may look like
1051 * \code
1052 * RA 12:31:56.17643
1053 * DEC 43:21:35.2531
1054 * F0 100.78634 1 0.00005
1055 * F1 2.34e-15
1056 * PEPOCH 54323.785634
1057 * \endcode
1058 * which shows several parameters mostly defined by the parameter name and a parameter value. However, the \c F0 value
1059 * contains 4 items. If a parameter has a \c 1 as the third entry then it means that this was a parameter that was fit
1060 * by TEMPO with the entry after the \c 1 being the 1 standard deviation error on that parameter. For parameters where
1061 * an error is present the code will attempt to search over that parameter using a Gaussian prior defined by the
1062 * 1 \f$ \sigma \f$ error value. Other parameters will be set as fixed by default. These can be overridden by the prior
1063 * file values described in \c initialise_prior().
1064 *
1065 * Based on the defined sky position defined in the par file a lookup table of the detector antenna response over time
1066 * and polarisation will be set by \c setup_lookup_tables().
1067 *
1068 * The function \c add_initial_variables() is used to pass the parameter values from the .par file to the algorithm.
1069 *
1070 * Using the parameters from the .par file the phase template, including the solar system and binary system barycentring
1071 * time delays will be setup. These define the phase template used to perform the initial heterodyne, which is used as
1072 * the reference in cases when phase parameters (other than the initial phase) are being searched over.
1073 *
1074 * Values used for scaling the parameters (to avoid dynamic range issues) are initialised although will be set as
1075 * default values.
1076 *
1077 * \param runState [in] A pointer to the LALInferenceRunState
1078 *
1079 * \sa setup_lookup_tables
1080 * \sa add_initial_variables
1081 * \sa get_phase_model
1082 * \sa add_correlation_matrix
1083 */
1085/* Read the PAR file of pulsar parameters and setup the code using them */
1086/* Generates lookup tables also */
1087{
1088 LALSource psr;
1089 PulsarParameters *pulsar = NULL;
1090 LALInferenceIFOData *data = runState->data;
1091 ProcessParamsTable *ppt = NULL;
1092 REAL8 DeltaT = 0.; /* maximum data time span */
1093
1094 ppt = LALInferenceGetProcParamVal( runState->commandLine, "--par-file" );
1095 if ( ppt == NULL ) {
1096 fprintf( stderr, "Must specify --par-file!\n" );
1097 exit( 1 );
1098 }
1099 CHAR *parFile = ppt->value;
1100
1101 /* get the pulsar parameters */
1102 pulsar = XLALReadTEMPOParFile( parFile );
1103
1104 REAL8 ra = 0.;
1105 if ( PulsarCheckParam( pulsar, "RA" ) ) {
1106 ra = PulsarGetREAL8Param( pulsar, "RA" );
1107 } else if ( PulsarCheckParam( pulsar, "RAJ" ) ) {
1108 ra = PulsarGetREAL8Param( pulsar, "RAJ" );
1109 } else {
1110 XLAL_ERROR_VOID( XLAL_EINVAL, "No source right ascension specified!" );
1111 }
1112 REAL8 dec = 0.;
1113 if ( PulsarCheckParam( pulsar, "DEC" ) ) {
1114 dec = PulsarGetREAL8Param( pulsar, "DEC" );
1115 } else if ( PulsarCheckParam( pulsar, "DECJ" ) ) {
1116 dec = PulsarGetREAL8Param( pulsar, "DECJ" );
1117 } else {
1118 XLAL_ERROR_VOID( XLAL_EINVAL, "No source declination specified!" );
1119 }
1120 psr.equatorialCoords.longitude = ra;
1121 psr.equatorialCoords.latitude = dec;
1123
1124 /* Setup lookup tables for amplitudes */
1125 setup_lookup_tables( runState, &psr );
1126
1127 runState->threads[0].model->params = XLALCalloc( 1, sizeof( LALInferenceVariables ) );
1128 runState->threads[0].model->domain = LAL_SIM_DOMAIN_TIME;
1129
1130 runState->threads[0].currentParams = XLALCalloc( 1, sizeof( LALInferenceVariables ) );
1131
1132 /* Add initial (unchanging) variables for the model from the par file */
1134
1135 /* check for binary model */
1136 CHAR *binarymodel = NULL;
1137 if ( LALInferenceCheckVariable( runState->threads[0].currentParams, "BINARY" ) ) {
1138 binarymodel = XLALStringDuplicate( *( CHAR ** )LALInferenceGetVariable( runState->threads[0].currentParams, "BINARY" ) );
1139
1140 /* now remove from runState->params (as it conflict with calls to LALInferenceCompareVariables in the proposal) */
1141 LALInferenceRemoveVariable( runState->threads[0].currentParams, "BINARY" );
1142 }
1143
1144 /* check for glitches */
1145 UINT4 glitches = 0;
1146 if ( LALInferenceCheckVariable( runState->threads[0].currentParams, "GLNUM" ) ) {
1147 glitches = 1;
1148 }
1149
1150 /* Setup barycentring delays */
1151 LALInferenceIFOModel *ifo_model = runState->threads[0].model->ifo;
1152 while ( data ) {
1153 REAL8Vector *freqFactors = NULL;
1154 UINT4 j = 0, k = 0;
1155 REAL8 dt = XLALGPSGetREAL8( &IFO_XTRA_DATA( ifo_model )->times->data[IFO_XTRA_DATA( ifo_model )->times->length - 1] ) - XLALGPSGetREAL8( &IFO_XTRA_DATA( ifo_model )->times->data[0] );
1156
1157 if ( dt > DeltaT ) {
1158 DeltaT = dt;
1159 }
1160
1161 freqFactors = *( REAL8Vector ** )LALInferenceGetVariable( ifo_model->params, "freqfactors" );
1162
1163 for ( j = 0; j < freqFactors->length; j++ ) {
1164 REAL8Vector *dts = NULL, *bdts = NULL, *glitchphase = NULL;
1165
1166 /* check whether using original Jones (2010) signal source model or a biaxial model (in the amplitude/phase parameterisation) */
1167 if ( freqFactors->length == 2 ) {
1168 UINT4 dummyvar = 1;
1169
1170 if ( LALInferenceGetProcParamVal( runState->commandLine, "--source-model" ) ) {
1171 LALInferenceAddVariable( ifo_model->params, "source_model", &dummyvar, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED );
1172 } else if ( LALInferenceGetProcParamVal( runState->commandLine, "--biaxial" ) ) {
1174 }
1175 }
1176
1177 /* add binary model to the general parameters */
1178 if ( binarymodel != NULL ) {
1179 LALInferenceAddVariable( ifo_model->params, "BINARY", &binarymodel, LALINFERENCE_string_t, LALINFERENCE_PARAM_FIXED );
1180 }
1181
1182 if ( glitches ) {
1183 LALInferenceAddVariable( ifo_model->params, "GLITCHES", &glitches, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED );
1184 }
1185
1186 if ( LALInferenceGetProcParamVal( runState->commandLine, "--inject-only" ) && LALInferenceGetProcParamVal( runState->commandLine, "--inject-coarse" ) ) {
1187 /* use this if just wanting to create an injection that has only been "coarse heterodyned" */
1188 dts = XLALCreateREAL8Vector( IFO_XTRA_DATA( ifo_model )->times->length );
1189 if ( binarymodel != NULL ) {
1190 bdts = XLALCreateREAL8Vector( IFO_XTRA_DATA( ifo_model )->times->length );
1191 }
1192 /* set the time delays to zero, so they do not get removed during the injected signal generation */
1193 for ( k = 0; k < dts->length; k++ ) {
1194 dts->data[k] = 0.;
1195 if ( binarymodel != NULL ) {
1196 bdts->data[k] = 0.;
1197 }
1198 }
1199 } else {
1200 dts = get_ssb_delay( pulsar, IFO_XTRA_DATA( ifo_model )->times, IFO_XTRA_DATA( ifo_model )->ephem, IFO_XTRA_DATA( ifo_model )->tdat, IFO_XTRA_DATA( ifo_model )->ttype, data->detector );
1201 bdts = get_bsb_delay( pulsar, IFO_XTRA_DATA( ifo_model )->times, dts, IFO_XTRA_DATA( ifo_model )->ephem );
1202 glitchphase = get_glitch_phase( pulsar, IFO_XTRA_DATA( ifo_model )->times, dts, bdts );
1203 }
1204
1206 if ( bdts != NULL ) {
1208 }
1209 if ( glitchphase != NULL ) {
1210 LALInferenceAddVariable( ifo_model->params, "glitch_phase", &glitchphase, LALINFERENCE_REAL8Vector_t, LALINFERENCE_PARAM_FIXED );
1211 }
1212
1213 data = data->next;
1214 ifo_model = ifo_model->next;
1215 }
1216 }
1217
1218 /* set frequency bin step from longest data time span */
1219 REAL8 df = 1. / ( 2.*DeltaT );
1221
1222 return;
1223}
1224
1225
1226/**
1227 * \brief Read in an ascii text file of nested samples, convert to posterior samples and create k-d tree
1228 *
1229 * This function reads in an ascii file of nested samples, converted them into posterior samples and them add them to a
1230 * k-d tree. The file name containing the samples must have been given as the command line argument \c sample-file and
1231 * there must be an accompanying file with the names of each column with the same file name with _params.txt appended.
1232 *
1233 * It is assumed that the samples are in ascending log likelihood order. It is also assumed that variable values in the
1234 * file (and are not likelihood-like values) are consistent with those given that have prior ranges defined in the prior
1235 * file/par file (as these ranges will be used as bounds in a k-d tree created from this data).
1236 *
1237 * As it is assumed that the points read in are from a previous nested sampling run the number of live points used for
1238 * that run are also required to be given with the \c sample-nlive argument. This will be used during the conversion to
1239 * posterior samples.
1240 *
1241 * If given the k-d tree cell size for using the posterior as a prior can be set with the \c prior-cell argument, if not
1242 * set this defaults to 32.
1243 *
1244 * In the future this will be altered so as to also read in an XML file of samples.
1245 *
1246 * NOTE: add the ability to read in multiple files and combine the posterior samples
1247 *
1248 * \param runState [in] A pointer to the LALInferenceRunState
1249 */
1251{
1252 ProcessParamsTable *ppt = NULL;
1253
1254 UINT4 Ncell = 8; /* default prior cell size */
1255
1256 UINT4 i = 0, k = 0, nsamps = 0, nnlive = 0, n = 0;
1257 UINT4Vector *nlive = NULL, *Nsamps = NULL;
1258 CHAR *nlivevals = NULL, *templives = NULL, *templive = NULL;
1259
1260 CHAR *sampfile = NULL;
1261 CHAR *tempsamps = NULL, *tempsamp = NULL;
1262
1263 LALStringVector *sampfilenames = NULL;
1264
1265 LALInferenceVariables ***params = NULL;
1266
1267 FILE *fp = NULL;
1268
1269 /* get names of nested sample file columns */
1270 ppt = LALInferenceGetProcParamVal( runState->commandLine, "--sample-files" );
1271 if ( ppt != NULL ) {
1272 sampfile = XLALStringDuplicate( ppt->value );
1273
1274 /* count the number of sample files from the comma separated vales and set their names */
1275 tempsamps = XLALStringDuplicate( sampfile );
1276
1277 nsamps = count_csv( tempsamps );
1278
1279 for ( i = 0; i < nsamps; i++ ) {
1280 tempsamp = XLALStringToken( &tempsamps, ",", 0 );
1281 sampfilenames = XLALAppendString2Vector( sampfilenames, tempsamp );
1282 }
1283 } else {
1284 return; /* no file so we don't use this function */
1285 }
1286
1287 ppt = LALInferenceGetProcParamVal( runState->commandLine, "--sample-nlives" );
1288 if ( ppt != NULL ) {
1289 nlivevals = XLALStringDuplicate( ppt->value );
1290
1291 templives = XLALStringDuplicate( nlivevals );
1292
1293 nnlive = count_csv( templives );
1294
1295 if ( nnlive != nsamps ) {
1296 XLAL_ERROR_VOID( XLAL_EINVAL, "Number of live points not equal to number of posterior files!" );
1297 }
1298
1299 for ( i = 0; i < nnlive; i++ ) {
1300 templive = XLALStringToken( &templives, ",", 0 );
1302 nlive->data[i] = atoi( templive );
1303 }
1304
1307 } else {
1308 fprintf( stderr, "Must set the number of live points used in the input nested samples file.\n\n" );
1309 fprintf( stderr, USAGE, runState->commandLine->program );
1310 exit( 0 );
1311 }
1312
1313 /* allocate memory for nested samples */
1314 params = XLALCalloc( nsamps, sizeof( LALInferenceVariables ** ) );
1315
1316 /* loop over files, convert to posterior samples and combine them */
1317 for ( n = 0; n < nsamps; n++ ) {
1318 CHAR *namefile = NULL, name[256];
1319 LALStringVector *paramNames = NULL;
1320
1321 i = 0;
1322
1323 /* initialise array as NULL */
1324 params[n] = NULL;
1325
1326 namefile = XLALStringDuplicate( sampfilenames->data[n] );
1327 namefile = XLALStringAppend( namefile, "_params.txt" );
1328
1329 /* check file exists */
1330 if ( fopen( namefile, "r" ) == NULL || fopen( sampfilenames->data[n], "r" ) == NULL ) {
1331 XLAL_ERROR_VOID( XLAL_EIO, "Cannot access either %s or %s!", namefile, sampfilenames->data[n] );
1332 }
1333
1334 /* read in parameter names */
1335 fp = fopen( namefile, "r" );
1336 while ( fscanf( fp, "%s", name ) != EOF ) {
1337 paramNames = XLALAppendString2Vector( paramNames, name );
1338 }
1339
1340 fclose( fp );
1341
1342 /* read in parameter values */
1343 fp = fopen( sampfilenames->data[n], "r" );
1344 while ( !feof( fp ) ) {
1345 REAL8 ps[paramNames->length];
1346 UINT4 j = 0;
1347
1348 for ( j = 0; j < paramNames->length; j++ ) {
1349 if ( fscanf( fp, "%lf", &ps[j] ) == EOF ) {
1350 break;
1351 }
1352 }
1353
1354 if ( feof( fp ) ) {
1355 break;
1356 }
1357
1358 /* dynamically allocate memory */
1359 params[n] = XLALRealloc( params[n], ( i + 1 ) * sizeof( LALInferenceVariables * ) );
1360 params[n][i] = NULL;
1361 params[n][i] = XLALCalloc( 1, sizeof( LALInferenceVariables ) );
1362
1363 /* add variables */
1364 for ( j = 0; j < paramNames->length; j++ ) {
1365 /* use vary type of this analyses parameters i.e. those set by the prior
1366 and par file, otherwise set the parameter to fixed */
1368
1369 if ( LALInferenceCheckVariable( runState->threads[0].currentParams, paramNames->data[j] ) ) {
1370 vary = LALInferenceGetVariableVaryType( runState->threads[0].currentParams, paramNames->data[j] );
1371 } else {
1373 }
1374
1375 LALInferenceAddVariable( params[n][i], paramNames->data[j], &ps[j], LALINFERENCE_REAL8_t, vary );
1376 }
1377
1378 i++;
1379 }
1380
1381 /* check that non-fixed, or output parameters actually do vary, otherwise
1382 complain */
1383 LALInferenceVariableItem *item1 = params[n][0]->head;
1384
1385 while ( item1 ) {
1386 UINT4 allsame = 0;
1387
1388 for ( k = 1; k < i; k++ ) {
1390
1391 if ( item1->vary != LALINFERENCE_PARAM_FIXED && item1->vary != LALINFERENCE_PARAM_OUTPUT ) {
1392 if ( *( REAL8 * )item1->value != *( REAL8 * )item2->value ) {
1393 allsame++;
1394 }
1395 }
1396 }
1397
1398 if ( ( item1->vary != LALINFERENCE_PARAM_FIXED && item1->vary != LALINFERENCE_PARAM_OUTPUT ) && allsame == 0 ) {
1399 XLAL_ERROR_VOID( XLAL_EFUNC, "Apparently variable parameter %s does not vary!\n", item1->name );
1400 }
1401
1402 item1 = item1->next;
1403 }
1404
1405 Nsamps = XLALResizeUINT4Vector( Nsamps, n + 1 );
1406 Nsamps->data[n] = i;
1407 }
1408
1413
1414 /* get cell size */
1415 ppt = LALInferenceGetProcParamVal( runState->commandLine, "--prior-cell" );
1416 if ( ppt != NULL ) {
1417 Ncell = atoi( ppt->value );
1418 }
1419
1420 LALInferenceAddVariable( runState->priorArgs, "kDTreePriorNcell", &Ncell, LALINFERENCE_UINT4_t,
1422
1423 /* convert samples to posterior */
1424 ns_to_posterior( runState );
1425
1426 /* create k-d tree of the samples for use as a prior */
1427 create_kdtree_prior( runState );
1428}
#define __func__
log an I/O error, i.e.
ProcessParamsTable * ppt
int j
int k
int PulsarCheckParam(const PulsarParameters *pars, const CHAR *name)
Check for the existence of the parameter name in the PulsarParameters structure.
PulsarParameters * XLALReadTEMPOParFile(const CHAR *pulsarAndPath)
Read in the parameters from a TEMPO(2) parameter file into a PulsarParameters structure.
REAL8 PulsarGetREAL8Param(const PulsarParameters *pars, const CHAR *name)
Return a REAL8 parameter.
const CHAR * PulsarGetStringParam(const PulsarParameters *pars, const CHAR *name)
Return a string parameter.
void PulsarFreeParams(PulsarParameters *pars)
Function to free memory from pulsar parameters.
REAL8 PulsarGetREAL8VectorParamIndividual(const PulsarParameters *pars, const CHAR *name)
Return an individual REAL8 value from a REAL8Vector parameter.
const char * name
Definition: SearchTiming.c:93
#define fscanf
#define fprintf
void XLALDestroyTokenList(TokenList *list)
int XLALCreateTokenList(TokenList **list, const CHAR *string, const CHAR *delimiters)
double e
char * XLALFileLoad(const char *path)
TimeCorrectionData * XLALInitTimeCorrections(const CHAR *timeCorrectionFile)
An XLAL interface for reading a time correction file containing a table of values for converting betw...
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
TimeCorrectionType
Enumerated type denoting the time system type to be produced in the solar system barycentring routine...
Definition: LALBarycenter.h:72
@ TIMECORRECTION_ORIGINAL
Definition: LALBarycenter.h:78
@ TIMECORRECTION_TCB
Definition: LALBarycenter.h:75
@ TIMECORRECTION_TDB
Definition: LALBarycenter.h:74
#define LAL_PI
double complex COMPLEX16
double REAL8
char CHAR
uint32_t UINT4
int32_t INT4
void LALInferenceAddVariable(LALInferenceVariables *vars, const char *name, const void *value, LALInferenceVariableType type, LALInferenceParamVaryType vary)
void LALInferenceRemoveVariable(LALInferenceVariables *vars, const char *name)
LALInferenceParamVaryType
void * LALInferenceGetVariable(const LALInferenceVariables *vars, const char *name)
ProcessParamsTable * LALInferenceGetProcParamVal(ProcessParamsTable *procparams, const char *name)
LALInferenceVariableItem * LALInferenceGetItem(const LALInferenceVariables *vars, const char *name)
#define DETNAMELEN
int LALInferenceCheckVariable(LALInferenceVariables *vars, const char *name)
LALInferenceParamVaryType LALInferenceGetVariableVaryType(LALInferenceVariables *vars, const char *name)
void LALInferenceSetVariable(LALInferenceVariables *vars, const char *name, const void *value)
LALINFERENCE_PARAM_OUTPUT
LALINFERENCE_PARAM_FIXED
LALINFERENCE_string_t
LALINFERENCE_REAL8_t
LALINFERENCE_INT4_t
LALINFERENCE_REAL8Vector_t
LALINFERENCE_void_ptr_t
LALINFERENCE_UINT4_t
LALINFERENCE_UINT4Vector_t
void * XLALMalloc(size_t n)
void * XLALCalloc(size_t m, size_t n)
void * XLALRealloc(void *p, size_t n)
void XLALFree(void *p)
LAL_SIM_DOMAIN_TIME
double XLALSimNoisePSDGEO(double f)
double XLALSimNoisePSDTAMA(double f)
double XLALSimNoisePSDaLIGOZeroDetHighPower(double f)
double XLALSimNoisePSDAdvVirgo(double f)
double XLALSimNoisePSDiLIGOSRD(double f)
double XLALSimNoisePSDVirgo(double f)
char char * XLALStringDuplicate(const char *s)
size_t XLALStringCopy(char *dst, const char *src, size_t size)
int char * XLALStringAppend(char *s, const char *append)
char * XLALStringToken(char **s, const char *delim, int empty)
RandomParams * XLALCreateRandomParams(INT4 seed)
void XLALDestroyRandomParams(RandomParams *params)
int XLALNormalDeviates(REAL4Vector *deviates, RandomParams *params)
const LALDetector * XLALGetSiteInfo(const CHAR *name)
Find the site geometry-information 'LALDetector' for given a detector name (or prefix).
Definition: SFTnaming.c:218
LIGOTimeGPSVector * XLALCreateTimestampVector(UINT4 len)
Allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:47
COORDINATESYSTEM_EQUATORIAL
LALStringVector * XLALAppendString2Vector(LALStringVector *vect, const CHAR *string)
COMPLEX16TimeSeries * XLALCreateCOMPLEX16TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
REAL8TimeSeries * XLALResizeREAL8TimeSeries(REAL8TimeSeries *series, int first, size_t length)
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
COMPLEX16TimeSeries * XLALResizeCOMPLEX16TimeSeries(COMPLEX16TimeSeries *series, int first, size_t length)
const LALUnit lalSecondUnit
REAL8Vector * XLALResizeREAL8Vector(REAL8Vector *vector, UINT4 length)
REAL4Vector * XLALCreateREAL4Vector(UINT4 length)
UINT4Vector * XLALResizeUINT4Vector(UINT4Vector *vector, UINT4 length)
void XLALDestroyREAL4Vector(REAL4Vector *vector)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
UINT4Vector * XLALCreateUINT4Vector(UINT4 length)
#define XLAL_ERROR_VOID(...)
#define XLAL_CHECK_VOID(assertion,...)
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EIO
XLAL_EINVAL
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
long long count
Definition: hwinject.c:371
float data[BLOCKSIZE]
Definition: hwinject.c:360
int gpstime
Definition: hwinject.c:365
n
tfile
nlive
void add_initial_variables(LALInferenceVariables *ini, PulsarParameters *pars)
Set up all the allowed variables for a known pulsar search This functions sets up all possible variab...
Definition: ppe_init.c:386
void setup_lookup_tables(LALInferenceRunState *runState, LALSource *source)
Sets the time angle antenna response lookup table.
Definition: ppe_init.c:296
void create_kdtree_prior(LALInferenceRunState *runState)
Create a k-d tree from prior samples.
void ns_to_posterior(LALInferenceRunState *runState)
Convert an array of nested samples to posterior samples.
REAL8Vector * get_glitch_phase(PulsarParameters *pars, LIGOTimeGPSVector *datatimes, REAL8Vector *dts, REAL8Vector *bdts)
Computes the phase from the glitch model.
Definition: ppe_models.c:798
REAL8Vector * get_bsb_delay(PulsarParameters *pars, LIGOTimeGPSVector *datatimes, REAL8Vector *dts, EphemerisData *edat)
Computes the delay between a pulsar in a binary system and the barycentre of the system.
Definition: ppe_models.c:762
REAL8Vector * get_ssb_delay(PulsarParameters *pars, LIGOTimeGPSVector *datatimes, EphemerisData *ephem, TimeCorrectionData *tdat, TimeCorrectionType ttype, LALDetector *detector)
Computes the delay between a GPS time at Earth and the solar system barycentre.
Definition: ppe_models.c:657
#define IFO_XTRA_DATA(ifo)
Definition: ppe_models.h:35
void samples_prior(LALInferenceRunState *runState)
Read in an ascii text file of nested samples, convert to posterior samples and create k-d tree.
void setup_from_par_file(LALInferenceRunState *runState)
Reads in the parameters of the pulsar being searched for.
void read_pulsar_data(LALInferenceRunState *runState)
Reads in the input gravitational wave data files, or creates fake data sets.
Definition: ppe_readdata.c:100
Header file for the data reading functions for the parameter estimation code for known pulsar searche...
void check_and_add_fixed_variable(LALInferenceVariables *vars, const char *name, void *value, LALInferenceVariableType type)
Add a variable, checking first if it already exists and is of type LALINFERENCE_PARAM_FIXED and if so...
Definition: ppe_utils.c:778
void compute_variance(LALInferenceIFOData *data, LALInferenceIFOModel *model)
Compute the noise variance for each data segment.
Definition: ppe_utils.c:38
UINT4Vector * get_chunk_lengths(LALInferenceIFOModel *ifo, UINT4 chunkMax)
Split the data into segments.
Definition: ppe_utils.c:98
TimeCorrectionType XLALAutoSetEphemerisFiles(CHAR **efile, CHAR **sfile, CHAR **tfile, PulsarParameters *pulsar, INT4 gpsstart, INT4 gpsend)
Automatically set the solar system ephemeris file based on environment variables and data time span.
Definition: ppe_utils.c:717
UINT4Vector * chop_n_merge(LALInferenceIFOData *data, UINT4 chunkMin, UINT4 chunkMax, UINT4 outputchunks)
Chops and remerges data into stationary segments.
Definition: ppe_utils.c:177
INT4 count_csv(CHAR *csvline)
Counts the number of comma separated values in a string.
Definition: ppe_utils.c:674
#define SQUARE(x)
#define MAXDETS
The maximum number of different detectors allowable.
#define FACTOR(x, y)
Macro that gives the integer number of times that x goes in to y.
#define CHUNKMAX
Default value of the maximum length into which the data can be split.
#define CHUNKMIN
Default value of the minimum length into which the data can be split.
LALDetector detectors[LAL_NUM_DETECTORS]
COMPLEX16Sequence * data
COMPLEX16 * data
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
REAL8TimeSeries * varTimeData
LALDetector * detector
char name[DETNAMELEN]
struct tagLALInferenceIFOData * next
COMPLEX16TimeSeries * compTimeData
struct tagLALInferenceIFOModel * next
LALDetector * detector
LALInferenceVariables * params
COMPLEX16TimeSeries * compTimeSignal
REAL8 * ifo_loglikelihoods
LALInferenceVariables * params
LALInferenceIFOModel * ifo
LALSimulationDomain domain
ProcessParamsTable * commandLine
LALInferenceVariables * algorithmParams
struct tagLALInferenceIFOData * data
LALInferenceVariables * priorArgs
LALInferenceThreadState * threads
LALInferenceVariables * currentParams
LALInferenceModel * model
char name[VARNAME_MAX]
struct tagVariableItem * next
LALInferenceParamVaryType vary
SkyPosition equatorialCoords
The PulsarParameters structure to contain a set of pulsar parameters.
REAL4 * data
REAL8Sequence * data
LIGOTimeGPS epoch
REAL8 * data
REAL8 longitude
REAL8 latitude
CoordinateSystem system
CHAR ** tokens
UINT4 nTokens
UINT4 * data
double df