Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInference 4.1.9.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALInferenceMultiBandTest.c
Go to the documentation of this file.
1/*
2 * LALInferenceBench.c: Benchmark LALInference functions
3 *
4 * Copyright (C) 2015 John Veitch
5 *
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with with program; see the file COPYING. If not, write to the
19 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
20 * MA 02110-1301 USA
21 */
22
23#include <stdio.h>
24#include <lal/LALInference.h>
25#include <lal/LALInferenceInit.h>
26#include <lal/LALInferenceReadData.h>
27#include <lal/LALInferenceCalibrationErrors.h>
28#include <lal/LALInferenceTemplate.h>
29#include <lal/LALInferenceLikelihood.h>
30#include <sys/resource.h>
31
32#ifdef __GNUC__
33#define UNUSED __attribute__ ((unused))
34#else
35#define UNUSED
36#endif
37
38const char HELPSTR[]=\
39"LALInferenceMultiBandTest: Unit test for consistency between multiband and regular template functions.\n\
40 Example (for 1.4-1.4 binary with seglen 32, srate 4096): \n\
41 $ ./LALInferenceMultiBandTest --psdlength 1000 --psdstart 1 --seglen 32 --srate 4096 --trigtime 0 --ifo H1 --H1-channel LALSimAdLIGO --H1-cache LALSimAdLIGO --dataseed 1324 --fix-chirpmass 1.218 --fix-q 1.0 --margphi\n\n\n\
42";
43
45
46
49{
50 return;
51}
52
55{
56 REAL8 oldPhase=0.0,mbPhase=0.0;
57 REAL8 h0MatchPlus,h0MatchCross;
58 REAL8 hmbMatchPlus,hmbMatchCross;
59
60 REAL8 target_snr = 50; /* Max SNR that we expect to handle */
61 REAL8 tolerance = 0.1/(target_snr*target_snr); /* Error in likelihood of 0.1 at target SNR */
62
63 COMPLEX16FrequencySeries *oldTemplatePlus=NULL,*oldTemplateCross=NULL,*newTemplatePlus=NULL,*newTemplateCross=NULL;
64
65
68 /* FIXME: Have to call the function once before result is repeatable!!! */
69 runState->likelihood(runState->threads[0].model->params,runState->data, runState->threads[0].model);
70 oldTemplatePlus = runState->threads[0].model->freqhPlus;
71 oldTemplateCross = runState->threads[0].model->freqhCross;
72
73 runState->threads[0].model->freqhPlus=XLALCreateCOMPLEX16FrequencySeries("mbtemplate",&oldTemplatePlus->epoch,oldTemplatePlus->f0,oldTemplatePlus->deltaF,&lalDimensionlessUnit,oldTemplatePlus->data->length);
74 runState->threads[0].model->freqhCross=XLALCreateCOMPLEX16FrequencySeries("mbtemplate",&oldTemplateCross->epoch,oldTemplateCross->f0,oldTemplateCross->deltaF,&lalDimensionlessUnit,oldTemplateCross->data->length);
75
76 /* Clear the template */
78
79 runState->likelihood(runState->threads[0].model->params,runState->data, runState->threads[0].model);
80 LALInferenceDumpWaveforms(runState->threads[0].model, "normal");
81 if(LALInferenceCheckVariable(runState->threads[0].model->params,"phase_maxl"))
82 oldPhase=LALInferenceGetREAL8Variable(runState->threads[0].model->params,"phase_maxl");
83
85
86 runState->likelihood(runState->threads[0].model->params,runState->data, runState->threads[0].model);
87 LALInferenceDumpWaveforms(runState->threads[0].model, "multiband");
88
89 REAL8 SNRsqPlus =LALInferenceComputeFrequencyDomainOverlap(runState->data,runState->threads[0].model->freqhPlus->data,runState->threads[0].model->freqhPlus->data);
90 REAL8 SNRsqCross =LALInferenceComputeFrequencyDomainOverlap(runState->data,runState->threads[0].model->freqhCross->data,runState->threads[0].model->freqhCross->data);
91
92 if(LALInferenceCheckVariable(runState->threads[0].model->params,"phase_maxl"))
93 mbPhase=LALInferenceGetREAL8Variable(runState->threads[0].model->params,"phase_maxl");
94
95 newTemplatePlus = runState->threads[0].model->freqhPlus;
96 newTemplateCross = runState->threads[0].model->freqhCross;
97
98 h0MatchPlus = LALInferenceComputeFrequencyDomainOverlap(runState->data,oldTemplatePlus->data,oldTemplatePlus->data);
99 h0MatchCross = LALInferenceComputeFrequencyDomainOverlap(runState->data,oldTemplateCross->data,oldTemplateCross->data);
100
101 hmbMatchPlus = LALInferenceComputeFrequencyDomainOverlap(runState->data,newTemplatePlus->data,newTemplatePlus->data);
102 hmbMatchCross = LALInferenceComputeFrequencyDomainOverlap(runState->data,newTemplateCross->data,newTemplateCross->data);
103 /* Want to check that <h0|h0>+<h_mb|h_mb>-2<h0|h_mb> < tolerance */
104
105 fprintf(stdout,"Parameter values:\n");
107
108 COMPLEX16 mismatchplus = compute_mismatch(runState->data, newTemplatePlus , oldTemplatePlus);
109 COMPLEX16 mismatchcross = compute_mismatch(runState->data, newTemplateCross , oldTemplateCross);
110
111 COMPLEX16 innerPlus = LALInferenceComputeFrequencyDomainComplexOverlap(runState->data, newTemplatePlus->data, oldTemplatePlus->data);
112 COMPLEX16 innerCross = LALInferenceComputeFrequencyDomainComplexOverlap(runState->data, newTemplateCross->data, oldTemplateCross->data);
113
114 int result = (1.0-cabs(innerPlus)/h0MatchPlus < tolerance) && (1.0 - cabs(innerCross)/h0MatchCross < tolerance);
115
116 fprintf(stdout,"\n\n");
117 fprintf(stdout,"|<h0|hmb>/<h0|h0>| = %lf (+) %lf (x)\n",(cabs(innerPlus)/h0MatchPlus),(cabs(innerCross)/h0MatchCross));
118 fprintf(stdout,"SNR = %lf (+), %lf (x)\n",sqrt(h0MatchPlus),sqrt(h0MatchCross));
119 fprintf(stdout,"SNR mb = %lf (+), %lf (x)\n",sqrt(hmbMatchPlus),sqrt(hmbMatchCross));
120 fprintf(stdout,"cplx <h|h0> (+): %lf*exp(%lfi), (x): %lf*exp(%lfi)\n",cabs(innerPlus),carg(innerPlus),cabs(innerCross),carg(innerCross));
121
122 fprintf(stdout,"Normalised mismatch (|h0-h|/|h0|)^2: %le (+), %le (x)\n",cabs(mismatchplus)/SNRsqPlus,cabs(mismatchcross)/SNRsqCross);
123 fprintf(stdout,"Tolerance = %le\n",tolerance);
124 fprintf(stdout,"Test result: %s\n",result?"passed":"failed");
125 fprintf(stdout,"Interpolation good up to SNR %lf (+), %lf (x)\n",
126 sqrt(tolerance*target_snr*target_snr *SNRsqPlus/ mismatchplus),
127 sqrt(tolerance*target_snr*target_snr *SNRsqCross/ mismatchplus) );
128 if(LALInferenceCheckVariable(runState->threads[0].model->params,"phase_maxl"))
129 fprintf(stdout,"max Like phase:\tnormal = %lf, phaseinterp = %lf\n",oldPhase,mbPhase);
130 return(result);
131}
132
133/* Computes <a-b|a-b> */
135{
140 return aa+bb-ab-ba;
141}
142
143int main(int argc, char *argv[]){
144
145 //#define default_command_line_len 18
146 //char * const default_command_line[default_command_line_len];
147
148 //"--psdlength","1000","--psdstart","1","--seglen","64","--srate","4096","--trigtime","0","--ifo","H1","--H1-channel","LALSimAdLIGO","--H1-cache","LALSimAdLIGO","--dataseed","1234"};
149
150
151 ProcessParamsTable *procParams = NULL;
152 LALInferenceRunState *runState=NULL;
153
154 procParams=LALInferenceParseCommandLine(argc,argv);
155
156 runState = LALInferenceInitRunState(procParams);
157
158 if(runState) {
160
161 /* Simulate calibration errors */
163 }
164
165 /* Set up the template and likelihood functions */
166 LALInferenceInitCBCThreads(runState,1);
168
169 /* Disable waveform caching */
170 runState->threads[0].model->waveformCache=NULL;
171
172 int result = compare_template(runState);
173
174 return(result);
175}
ProcessParamsTable * LALInferenceParseCommandLine(int argc, char *argv[])
void LALInferenceApplyCalibrationErrors(LALInferenceIFOData *IFOdata, ProcessParamsTable *commandLine)
LALInferenceRunState * LALInferenceInitRunState(ProcessParamsTable *command_line)
void LALInferenceInitCBCThreads(LALInferenceRunState *run_state, INT4 nthreads)
int main(int argc, char *argv[])
void LALInferenceTemplateNoop(UNUSED LALInferenceModel *model)
COMPLEX16 compute_mismatch(LALInferenceIFOData *data, COMPLEX16FrequencySeries *a, COMPLEX16FrequencySeries *b)
int compare_template(LALInferenceRunState *runState)
const char HELPSTR[]
void LALInferenceInjectInspiralSignal(LALInferenceIFOData *IFOdata, ProcessParamsTable *commandLine)
#define fprintf
sigmaKerr data[0]
COMPLEX16FrequencySeries * XLALCreateCOMPLEX16FrequencySeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaF, const LALUnit *sampleUnits, size_t length)
double complex COMPLEX16
double REAL8
REAL8 LALInferenceGetREAL8Variable(LALInferenceVariables *vars, const char *name)
void LALInferenceDumpWaveforms(LALInferenceModel *model, const char *basefilename)
Dump all waveforms from the ifodata structure.
void LALInferencePrintVariables(LALInferenceVariables *var)
output contents of a 'LALInferenceVariables' structure * / / * (by now only prints names and types,...
Definition: LALInference.c:798
int LALInferenceCheckVariable(LALInferenceVariables *vars, const char *name)
Checks for name being present in vars returns 1(==true) or 0.
Definition: LALInference.c:525
void LALInferenceInitLikelihood(LALInferenceRunState *runState)
Initialisation function which reads runState->commaneLine and sets up the likelihood function accordi...
COMPLEX16 LALInferenceComputeFrequencyDomainComplexOverlap(LALInferenceIFOData *dataPtr, COMPLEX16Vector *freqData1, COMPLEX16Vector *freqData2)
Computes the complex <x|y> overlap.
REAL8 LALInferenceComputeFrequencyDomainOverlap(LALInferenceIFOData *dataPtr, COMPLEX16Vector *freqData1, COMPLEX16Vector *freqData2)
Computes the <x|y> overlap in the Fourier domain.
void LALInferenceTemplateNullFreqdomain(LALInferenceModel *model)
Returns a frequency-domain 'null' template (all zeroes, implying no signal present).
void LALInferenceTemplateXLALSimInspiralChooseWaveformPhaseInterpolated(LALInferenceModel *model)
void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
"XLALSimInspiralChooseWaveform{TD,FD}" wrapper.
static const INT4 a
const LALUnit lalDimensionlessUnit
COMPLEX16Sequence * data
Structure to contain IFO data.
Definition: LALInference.h:625
Structure to constain a model and its parameters.
Definition: LALInference.h:436
COMPLEX16FrequencySeries * freqhPlus
Time series model buffers.
Definition: LALInference.h:454
LALSimInspiralWaveformCache * waveformCache
Definition: LALInference.h:458
LALInferenceVariables * params
Definition: LALInference.h:437
COMPLEX16FrequencySeries * freqhCross
Definition: LALInference.h:454
LALInferenceTemplateFunction templt
Domain of model.
Definition: LALInference.h:440
Structure containing inference run state This includes pointers to the function types required to run...
Definition: LALInference.h:592
ProcessParamsTable * commandLine
Definition: LALInference.h:593
struct tagLALInferenceIFOData * data
Log sample, i.e.
Definition: LALInference.h:601
LALInferenceThreadState * threads
Definition: LALInference.h:612
LALInferenceLikelihoodFunction likelihood
MultiNest prior for the parameters.
Definition: LALInference.h:599
LALInferenceModel * model
Cycle of proposals to call.
Definition: LALInference.h:548