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
pulsar_parameter_estimation_nested.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2012 Matthew Pitkin, Colin Gill, John Veitch
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/* functions to create the likelihood for a pulsar search to be used with the
21LALInference tools */
22
23/**
24 * \file
25 * \ingroup lalpulsar_bin_HeterodyneSearch
26 * \author Matthew Pitkin, John Veitch, Max Isi, Colin Gill
27 *
28 * \brief Parameter estimation code for known pulsar searches using the nested sampling algorithm.
29 *
30 * ### Description ###
31 *
32 * This code is used to perform parameter estimation and evidence calculation in targeted/semi-targeted searches for
33 * gravitational waves from known pulsars. It may also be used to follow-up on signal candidates from semi-coherent all-sky
34 * searches for unknown sources.
35 *
36 * It uses the Bayesian technique of 'Nested Sampling' to sample over a defined prior parameter space (unknown signal
37 * parameters such as the gravitational wave amplitude). These samples can then be used to create posterior probability
38 * density functions of the unknown parameters. The samples are also used to calculate the Bayesian evidence, also known as
39 * the marginal likelihood, for a given signal model. This can be compared with other models, in particular the model that
40 * the data is described by Gaussian noise alone.
41 *
42 * As input the code requires time domain data that has been heterodyned using the known (or close to) phase evolution of
43 * the pulsar. The time domain input should consist of a three column text file containing the GPS time stamp of the data
44 * point, the real part of the heterodyned data and the imaginary part of the heterodyned data, e.g.
45 * \code
46 * 900000000.000000 1.867532e-24 -7.675329e-25
47 * 900000060.000000 2.783651e-24 3.654386e-25
48 * ...
49 * \endcode
50 *
51 * Most commonly such data will have a sample rate of 1/60 Hz, giving a bandwidth of the same amount, but the code can
52 * accept any rate.
53 *
54 * The code also requires that you specify which parameters are to be searched over, and the prior ranges over these. Any
55 * of the signal parameters can be searched over, including frequency, sky position and binary system parameters, although
56 * the bandwidth of the data and search efficiency need to be taken into account.
57 *
58 * The 'Nested Sampling' algorithm (developed by \cite Skilling2006) used is that defined in \c LALinferenceNestedSampler.c
59 * (see \cite VeitchVecchio2010). It is essentially an efficient way to perform the integral
60 * \f[
61 * Z = \int^{\mathbf{\theta}} p(d|\mathbf{\theta}) p(\mathbf{\theta}) \mathrm{d}\mathbf{\theta},
62 * \f]
63 * where \f$ \mathbf{\theta} \f$ is a vector of parameters, \f$ p(d|\mathbf{\theta}) \f$ is the likelihood of the data
64 * given the parameters, and \f$ p(\mathbf{\theta}) \f$ is the prior on the parameters. It does this by changing the
65 * multi-dimensional integral over N parameters into a one-dimensional integral
66 * \f[
67 * Z = \int^X L(X) \mathrm{d}X \approx \sum_i L(X_i) \Delta{}X_i,
68 * \f]
69 * where \f$ L(X) \f$ is the likelihood, and \f$ X \f$ is the prior mass. The algorithm will draw a number ( \f$ N \f$ ) of
70 * samples (live points) from the parameter priors, calculate the likelihood for each point and find the lowest likelihood
71 * value. The lowest likelihood value will be added to the summation in the above equation, with \f$ \log{\Delta{}X_i}
72 * \approx 1/N \f$ coming from the fact that the prior would be normalised to unity and therefore each point should occupy
73 * an equal fraction and at each iteration the prior volume will decrease geometrically (for \f$ \log{\Delta{}X_0} = 0 \f$ ).
74 * A new point is then drawn from the prior with the criterion that it has a higher likelihood than the previous lowest
75 * point and substitutes that point. To draw the new point a Markov Chain Monte Carlo (MCMC) procedure is used. The procedure
76 * is continued until a stopping criterion is reached, which in this case is that the remaining prior volume is less than the
77 * \c tolerance value set (see below). The implementation of this can be seen in \cite VeitchVecchio2010 .
78 *
79 * ### Usage ###
80 *
81 * The usage format is given below and can also be found by running the code with
82 * \code
83 * lalpulsar_parameter_estimation_nested --help
84 * \endcode
85 *
86 * An example of running the code to search over the four unknown parameters \f$ h_0 \f$ , \f$ \phi_0 \f$ , \f$ \psi \f$
87 * and \f$ \cos{\iota} \f$ , for pulsar J0534-2200, given heterodyned time domain data from the H1 detector in the file
88 * \c finehet_J0534-2200_H1, is:
89 * \code
90 * lalpulsar_parameter_estimation_nested --detectors H1 --par-file J0534-2200.par --input-files finehet_J0534-2200_H1 --outfile ns_J0534-2200.hdf --prior-file prior_J0534-2200.txt --ephem-earth lscsoft/share/lalpulsar/earth05-09.dat --ephem-sun lscsoft/share/lalpulsar/sun05-09.dat --Nlive 1000 --Nmcmcinitial 0 --tolerance 0.25
91 * \endcode
92 * The \c par-file is a TEMPO(2)-style file containing the parameters of the pulsar used to perform the heterodyne (the
93 * frequency parameters are the rotation frequency and therefore not necessarily the gravitational wave frequency) e.g.
94 * \code
95 * RA 12:54:31.87523895
96 * DEC -54:12:43.6572033
97 * PMRA 1.7
98 * PMDEC 2.8
99 * POSEPOCH 54320.8531
100 * F0 123.7896438753
101 * F1 4.592e-15
102 * PEPOCH 54324.8753
103 * \endcode
104 * The \c prior-file is a text file containing a list of the parameters to be searched over, the prior type ("uniform" or
105 * "gaussian") and their given lower/mean and upper/standard deviation ranges e.g.
106 * \code
107 * h0 uniform 0 1e-21
108 * phi0 uniform 0 6.283185307179586
109 * cosiota uniform -1 1
110 * psi uniform -0.785398163397448 0.785398163397448
111 * \endcode
112 * Note that if searching over frequency parameters the ranges specified in the \c prior-file should be given in terms of
113 * the pulsars rotation frequency and not necessarily the gravitational wave frequency e.g. for a triaxial star emitting
114 * gravitational waves at 100 Hz (which will be at twice the rotation frequency) if you wanted to search over 99.999 to
115 * 100.001 Hz then you should used
116 * \code
117 * f0 uniform 49.9995 50.0005
118 * \endcode
119 *
120 * An example of running the code as above, but this time on fake data created using the Advanced LIGO design noise curves
121 * and with a signal injected into the data is:
122 * \code
123 * lalpulsar_parameter_estimation_nested --fake-data AH1 --inject-file fake.par --par-file fake.par --outfile ns_fake.hdf --prior-file prior_fake.txt --ephem-earth lscsoft/share/lalpulsar/earth05-09.dat --ephem-sun lscsoft/share/lalpulsar/sun05-09.dat --Nlive 1000 --Nmcmcinitial 0 --tolerance 0.25
124 * \endcode
125 * In this case the \c inject-file parameter file must contain the values of \c h0, \c phi0, \c psi and \c cosiota,
126 * otherwise these will be set to zero by default. The parameter files given for \c inject-file and \c par-file do not
127 * have to be the same - the injection can be offset from the 'heterodyned' values, which will be reflected in the data. If
128 * an \c inject-output file is also specified then the fake data containing the signal, and a fake signal-only data set,
129 * will be output.
130 */
131
132#include "config.h"
134#include "ppe_utils.h"
135#include "ppe_init.h"
136#include "ppe_readdata.h"
137#include "ppe_inject.h"
138#include "ppe_models.h"
139#include "ppe_likelihood.h"
140#include "ppe_testing.h"
141#include "ppe_roq.h"
142
143/* global variables */
145
146INT4 main( INT4 argc, CHAR *argv[] )
147{
148 ProcessParamsTable *param_table, *testgausslike;
149 LALInferenceRunState runState;
150 REAL8 logZnoise = 0.;
151 struct timeval time1, time2;
152 gettimeofday( &time1, NULL ); /* time program */
153
154 /* set error handler to abort in main function */
156
157 /* Get ProcParamsTable from input arguments */
159 runState.commandLine = param_table;
160
161 /* Initialise the algorithm structures from the command line arguments */
162 /* Include setting up random number generator etc */
163 initialise_algorithm( &runState );
164
165 /* check if testing with hardcoded Gaussian likelihood */
166 testgausslike = LALInferenceGetProcParamVal( param_table, "--test-gaussian-likelihood" );
167
168 /* read in data */
169 if ( !testgausslike ) {
170 read_pulsar_data( &runState );
171 } else {
172 /* initialise some required values if running on test Gaussian likelihood */
173 REAL8 h0val = 0., h0sigma = 2.5e-24, h0mean = 0.;
174 runState.data = NULL;
175 runState.threads[0].model = XLALCalloc( 1, sizeof( LALInferenceModel ) );
176 runState.data = XLALCalloc( 1, sizeof( LALInferenceIFOData ) );
177 runState.data->likeli_counter = 0;
178 runState.data->templa_counter = 0;
179 runState.data->next = NULL;
180 runState.threads[0].model->ifo = XLALMalloc( sizeof( LALInferenceIFOModel ) );
181 runState.threads[0].model->ifo->extraData = XLALCalloc( 1, sizeof( IFOModelExtraData ) );
182 runState.threads[0].model->ifo->params = XLALCalloc( 1, sizeof( LALInferenceVariables ) );
183 runState.threads[0].model->ifo->next = NULL;
184 runState.threads[0].model->ifo_loglikelihoods = XLALMalloc( sizeof( REAL8 ) );
185 runState.threads[0].model->ifo_SNRs = XLALMalloc( sizeof( REAL8 ) );
186 runState.threads[0].currentParams = XLALCalloc( 1, sizeof( LALInferenceVariables ) );
187
188 /* add parameters of the Gaussian */
189 LALInferenceAddVariable( runState.threads[0].currentParams, "H0", &h0val, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED ); /* add H0 parameter */
190 if ( LALInferenceGetProcParamVal( param_table, "--test-gaussian-sigma" ) ) {
191 h0sigma = atof( LALInferenceGetProcParamVal( param_table, "--test-gaussian-sigma" )->value );
192 }
193 LALInferenceAddVariable( runState.threads[0].currentParams, "H0SIGMA", &h0sigma, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED ); /* add standard deviation */
194 if ( LALInferenceGetProcParamVal( param_table, "--test-gaussian-mean" ) ) {
195 h0mean = atof( LALInferenceGetProcParamVal( param_table, "--test-gaussian-mean" )->value );
196 }
197 LALInferenceAddVariable( runState.threads[0].currentParams, "H0MEAN", &h0mean, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED ); /* add standard deviation */
198 }
199
200 /* set algorithm to use Nested Sampling */
203
204 /* set likelihood function */
205 if ( testgausslike ) {
206 runState.likelihood = &test_gaussian_log_likelihood; /* run on a simple test Gaussian likelihood */
207 } else {
209 }
210
211 /* set prior function */
212 runState.prior = &priorFunction;
213
214 /* Generate the lookup tables and read parameters from par file */
215 if ( !testgausslike ) {
216 setup_from_par_file( &runState );
217 }
218
219 /* set signal model/template */
220 if ( !testgausslike ) {
221 runState.threads[0].model->templt = &get_pulsar_model;
222 }
223
224 /* add injections if requested */
225 if ( !testgausslike ) {
226 inject_signal( &runState );
227 }
228
229 /* Initialise the prior distribution given the command line arguments */
230 initialise_prior( &runState );
231
232 /* create sum square of the data to speed up the likelihood calculation */
233 if ( !testgausslike ) {
234 sum_data( &runState );
235 }
236
237 /* check whether using reduced order quadrature */
238 if ( !testgausslike ) {
239 generate_interpolant( &runState );
240 }
241
242 if ( !testgausslike ) {
243 gridOutput( &runState );
244 }
245
246 /* get noise likelihood and add as variable to runState */
247 if ( !testgausslike ) {
248 logZnoise = noise_only_likelihood( &runState );
249 } else {
250 logZnoise = 0.;
251 }
253
254 /* Create live points array and fill initial parameters */
255 if ( !LALInferenceGetProcParamVal( param_table, "--compare-likelihoods" ) ) {
257 }
258
259 /* output the live points sampled from the prior */
260 outputPriorSamples( &runState );
261
262 /* Initialise the MCMC proposal distribution */
263 initialise_proposal( &runState );
264
265 /* Set up threads */
266 initialise_threads( &runState, 1 );
267
268 if ( !LALInferenceGetProcParamVal( param_table, "--compare-likelihoods" ) ) {
269 /* Call the nested sampling algorithm */
270 runState.algorithm( &runState );
271 } else {
272 /* compare likelihoods from previous run */
273 compare_likelihoods( &runState );
274 return 0;
275 }
276
277 /* get SNR of highest likelihood point */
278 if ( !testgausslike ) {
279 get_loudest_snr( &runState );
280 }
281
282 /* output log evidence and 95% upper limit of test Gaussian likelihood */
283 if ( testgausslike ) {
284 test_gaussian_output( &runState );
285 }
286
287 /* close timing file */
288 if ( LALInferenceCheckVariable( runState.algorithmParams, "timefile" ) ) {
289 gettimeofday( &time2, NULL );
290 REAL8 tottime = ( REAL8 )( ( time2.tv_sec + time2.tv_usec * 1.e-6 ) - ( time1.tv_sec + time1.tv_usec * 1.e-6 ) );
291 FILE *timefile = *( FILE ** )LALInferenceGetVariable( runState.algorithmParams, "timefile" );
292 UINT4 timenum = *( UINT4 * )LALInferenceGetVariable( runState.algorithmParams, "timenum" );
293 fprintf( timefile, "[%d] %s: %.9le secs\n", timenum, __func__, tottime );
294 fclose( timefile );
295 }
296
297 return 0;
298}
#define __func__
log an I/O error, i.e.
#define fprintf
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 * LALInferenceGetVariable(const LALInferenceVariables *vars, const char *name)
ProcessParamsTable * LALInferenceGetProcParamVal(ProcessParamsTable *procparams, const char *name)
ProcessParamsTable * LALInferenceParseCommandLine(int argc, char **argv)
int LALInferenceCheckVariable(LALInferenceVariables *vars, const char *name)
LALINFERENCE_PARAM_FIXED
LALINFERENCE_REAL8_t
INT4 LALInferenceNestedSamplingOneStep(LALInferenceRunState *runState)
void * XLALMalloc(size_t n)
void * XLALCalloc(size_t m, size_t n)
XLALErrorHandlerType * XLALSetErrorHandler(XLALErrorHandlerType *newHandler)
void XLALExitErrorHandler(const char *func, const char *file, int line, int errnum)
list param_table
void initialise_proposal(LALInferenceRunState *runState)
Initialise the MCMC proposal distribution for sampling new points.
Definition: ppe_init.c:918
void initialise_prior(LALInferenceRunState *runState)
Sets up the parameters to be searched over and their prior ranges.
Definition: ppe_init.c:620
void sum_data(LALInferenceRunState *runState)
Calculates the sum of the square of the data and model terms.
Definition: ppe_init.c:1110
void initialise_threads(LALInferenceRunState *state, INT4 nthreads)
Definition: ppe_init.c:1719
void initialise_algorithm(LALInferenceRunState *runState)
Initialises the nested sampling algorithm control.
Definition: ppe_init.c:138
void setup_live_points_array_wrapper(LALInferenceRunState *runState)
A wrapper around LALInferenceSetupLivePointsArray.
Definition: ppe_init.c:100
void nested_sampling_algorithm_wrapper(LALInferenceRunState *runState)
A wrapper around LALInferenceNestedSamplingAlgorithm.
Definition: ppe_init.c:57
Header file for the initialisation functions for the parameter estimation code for known pulsar searc...
void inject_signal(LALInferenceRunState *runState)
Inject a simulated signal into the data.
Definition: ppe_inject.c:47
void get_loudest_snr(LALInferenceRunState *runState)
Get the signal-to-noise ratio of the maximum likelihood signal.
Definition: ppe_inject.c:647
Header file for the signal injection functions for the parameter estimation code for known pulsar sea...
REAL8 priorFunction(LALInferenceRunState *runState, LALInferenceVariables *params, UNUSED LALInferenceModel *model)
The prior function.
REAL8 noise_only_likelihood(LALInferenceRunState *runState)
Calculate the natural logarithm of the evidence that the data consists of only Gaussian noise The fun...
REAL8 pulsar_log_likelihood(LALInferenceVariables *vars, LALInferenceIFOData *data, LALInferenceModel *get_model)
The log likelihood function.
Header file for the likelihood and prior functions used in parameter estimation code for known pulsar...
void get_pulsar_model(LALInferenceModel *model)
Defines the pulsar model/template to use.
Definition: ppe_models.c:52
Header file for the signal models functions used in parameter estimation code for known pulsar search...
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 generate_interpolant(LALInferenceRunState *runState)
Generate an orthonormal basis set of waveforms from a training set.
Definition: ppe_roq.c:82
Header file for the reduced order quadrature generation used in parameter estimation code for known p...
void test_gaussian_output(LALInferenceRunState *runState)
Output the analytic evidence for the test Gaussian likelihood and a 95% upper limit.
Definition: ppe_testing.c:245
REAL8 test_gaussian_log_likelihood(LALInferenceVariables *vars, LALInferenceIFOData *data, LALInferenceModel *get_model)
Test the sampler using a Gaussian likelihood.
Definition: ppe_testing.c:209
void compare_likelihoods(LALInferenceRunState *rs)
Read in an ascii text file of nested samples and compare the log likelihoods.
Definition: ppe_testing.c:450
void outputPriorSamples(LALInferenceRunState *runState)
Output a number of prior samples based on the initial live points.
Definition: ppe_testing.c:385
void gridOutput(LALInferenceRunState *runState)
A test function to calculate a 1D posterior on a grid.
Definition: ppe_testing.c:51
Header file for functions used in testing the parameter estimation code for known pulsar searches usi...
Header file for the helper functions for the parameter estimation code for known pulsar searches usin...
LALStringVector * corlist
INT4 main(INT4 argc, CHAR *argv[])
Header file for the parameter estimation code for known pulsar searches using the nested sampling alg...
struct tagLALInferenceIFOModel * next
LALInferenceVariables * params
REAL8 * ifo_loglikelihoods
LALInferenceIFOModel * ifo
LALInferenceTemplateFunction templt
ProcessParamsTable * commandLine
LALInferenceVariables * algorithmParams
struct tagLALInferenceIFOData * data
LALInferenceThreadState * threads
LALInferenceAlgorithm algorithm
LALInferenceEvolveOneStepFunction evolve
LALInferencePriorFunction prior
LALInferenceLikelihoodFunction likelihood
LALInferenceVariables * currentParams
LALInferenceModel * model