LALInferenceInitCBC.c 114 KB
Newer Older
1
2
3
/*
 *  LALInferenceCBCInit.c:  Bayesian Followup initialisation routines.
 *
4
 *  Copyright (C) 2012 Vivien Raymond, John Veitch, Salvatore Vitale
5
6
7
8
9
10
11
12
13
14
15
16
17
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with with program; see the file COPYING. If not, write to the
18
19
 *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 *  MA  02110-1301  USA
20
21
22
23
 */


#include <stdio.h>
24
#include <assert.h>
25
#include <errno.h>
26
27
28
29
30
31
32
33
34
35
36
37
38
#include <lal/Date.h>
#include <lal/GenerateInspiral.h>
#include <lal/LALInference.h>
#include <lal/FrequencySeries.h>
#include <lal/Units.h>
#include <lal/StringInput.h>
#include <lal/LIGOLwXMLInspiralRead.h>
#include <lal/TimeSeries.h>
#include <lal/LALInferencePrior.h>
#include <lal/LALInferenceTemplate.h>
#include <lal/LALInferenceProposal.h>
#include <lal/LALInferenceLikelihood.h>
#include <lal/LALInferenceReadData.h>
39
#include <lal/LALInferenceInit.h>
40
#include <lal/LALInferenceCalibrationErrors.h>
John Douglas Veitch's avatar
John Douglas Veitch committed
41
#include <lal/LALSimNeutronStar.h>
42

43
44
45
46
47
48
49
static int checkParamInList(const char *list, const char *param);
static int checkParamInList(const char *list, const char *param)
{
  /* Check for param in comma-seperated list */
  char *post=NULL,*pos=NULL;
  if (list==NULL) return 0;
  if (param==NULL) return 0;
50

51
  if(!(pos=strstr(list,param))) return 0;
52

53
54
55
56
57
  /* The string is a substring. Check that it is a token */
  /* Check the character before and after */
  if(pos!=list)
  if(*(pos-1)!=',')
  return 0;
58

59
60
61
62
63
64
  post=&(pos[strlen(param)]);
  if(*post!='\0')
  if(*post!=',')
  return 0;
  return 1;
}
65
66

static void print_flags_orders_warning(SimInspiralTable *injt, ProcessParamsTable *commline);
67
static void LALInferenceInitSpinVariables(LALInferenceRunState *state, LALInferenceModel *model);
68
static void LALInferenceInitMassVariables(LALInferenceRunState *state);
69
static void LALInferenceInitNonGRParams(LALInferenceRunState *state, LALInferenceModel *model);
70
static void LALInferenceCheckApproximantNeeds(LALInferenceRunState *state,Approximant approx);
71

72
73
74

/* Initialize a bare-bones run-state,
   calling the "ReadData()" function to gather data & PSD from files,
75
   sets up the random seed and rng, and initializes other variables accordingly.
76
*/
77
LALInferenceRunState *LALInferenceInitRunState(ProcessParamsTable *command_line) {
78
79
80
81
    ProcessParamsTable *ppt=NULL;
    INT4 randomseed;
    FILE *devrandom;
    struct timeval tv;
82
83
84
85
86
87
    LALInferenceRunState *run_state = XLALCalloc(1, sizeof(LALInferenceRunState));

    /* Check that command line is consistent first */
    LALInferenceCheckOptionsConsistency(command_line);
    run_state->commandLine = command_line;

88
89
90
91
92
    /* Initialize parameters structure */
    run_state->algorithmParams = XLALCalloc(1, sizeof(LALInferenceVariables));
    run_state->priorArgs = XLALCalloc(1, sizeof(LALInferenceVariables));
    run_state->proposalArgs = XLALCalloc(1, sizeof(LALInferenceVariables));

93
94
95
96
97
    /* Read data from files or generate fake data */
    run_state->data = LALInferenceReadData(command_line);
    if (run_state->data == NULL)
        return(NULL);

98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
    /* Setup the random number generator */
    gsl_rng_env_setup();
    run_state->GSLrandom = gsl_rng_alloc(gsl_rng_mt19937);

    /* Use clocktime if seed isn't provided */
    ppt = LALInferenceGetProcParamVal(command_line, "--randomseed");
    if (ppt)
        randomseed = atoi(ppt->value);
    else {
        if ((devrandom = fopen("/dev/urandom","r")) == NULL) {
            gettimeofday(&tv, 0);
            randomseed = tv.tv_sec + tv.tv_usec;
        } else {
            fread(&randomseed, sizeof(randomseed), 1, devrandom);
            fclose(devrandom);
        }
    }
    gsl_rng_set(run_state->GSLrandom, randomseed);

    /* Save the random seed */
    LALInferenceAddVariable(run_state->algorithmParams, "random_seed", &randomseed,
                            LALINFERENCE_INT4_t, LALINFERENCE_PARAM_OUTPUT);
120
121
122
123

    return(run_state);
}

124
125
/* Draw initial parameters for each of the threads in run state */
void LALInferenceDrawThreads(LALInferenceRunState *run_state) {
Vivien Raymond's avatar
Vivien Raymond committed
126
127
    if (run_state == NULL)
        return;
Ben Farr's avatar
Ben Farr committed
128
129
130

    if (run_state->threads == NULL) {
        XLALPrintError("Error: LALInferenceDrawThreads expects initialized run_state->threads\n");
131
        XLAL_ERROR_VOID(XLAL_EINVAL);
Ben Farr's avatar
Ben Farr committed
132
133
    }

134
    LALInferenceThreadState *thread = &(run_state->threads[0]);
135
136
137
138
139
    INT4 t;

    /* If using a malmquist prior, force a strict prior window on distance for starting point, otherwise
     * the approximate prior draws are very unlikely to be within the malmquist prior */
    REAL8 dist_low, dist_high;
Ben Farr's avatar
Ben Farr committed
140
141
    REAL8 restricted_dist_low = log(10.0);
    REAL8 restricted_dist_high = log(100.0);
142
143
    INT4 changed_dist = 0;
    if (LALInferenceCheckVariable(run_state->priorArgs, "malmquist") &&
Ben Farr's avatar
Ben Farr committed
144
        LALInferenceCheckVariableNonFixed(thread->currentParams, "logdistance")) {
145
        changed_dist = 1;
Ben Farr's avatar
Ben Farr committed
146
147
148
        LALInferenceGetMinMaxPrior(run_state->priorArgs, "logdistance", &dist_low, &dist_high);
        LALInferenceRemoveMinMaxPrior(run_state->priorArgs, "logdistance");
        LALInferenceAddMinMaxPrior(run_state->priorArgs, "logdistance",
149
150
151
152
153
154
155
                                   &restricted_dist_low, &restricted_dist_high, LALINFERENCE_REAL8_t);
    }

    /* If the currentParams are not in the prior, overwrite and pick paramaters
     *   from the priors. OVERWRITE EVEN USER CHOICES.
     *   (necessary for complicated prior shapes where
     *   LALInferenceCyclicReflectiveBound() is not enough) */
156
    #pragma omp parallel for private(thread)
157
    for (t = 0; t < run_state->nthreads; t++) {
Ben Farr's avatar
Ben Farr committed
158
159
        LALInferenceVariables *priorDraw = XLALCalloc(1, sizeof(LALInferenceVariables));

160
        thread = &(run_state->threads[t]);
Vivien Raymond's avatar
Vivien Raymond committed
161

162
163
164
165
166
        /* Try not to clobber values given on the command line */
        LALInferenceCopyVariables(thread->currentParams, priorDraw);
        LALInferenceDrawApproxPrior(thread, priorDraw, priorDraw);
        LALInferenceCopyUnsetREAL8Variables(priorDraw, thread->currentParams,
                                            run_state->commandLine);
Vivien Raymond's avatar
Vivien Raymond committed
167

Ben Farr's avatar
Ben Farr committed
168
169
170
        while(isinf(run_state->prior(run_state,
                                     thread->currentParams,
                                     thread->model))) {
171
172
173
174
175
176
177
178
179
180
            LALInferenceDrawApproxPrior(thread,
                                        thread->currentParams,
                                        thread->currentParams);
        }

        /* Make sure that our initial value is within the
        *     prior-supported volume. */
        LALInferenceCyclicReflectiveBound(thread->currentParams, run_state->priorArgs);

        /* Initialize starting likelihood and prior */
Ben Farr's avatar
Ben Farr committed
181
182
183
        thread->currentPrior = run_state->prior(run_state,
                                                thread->currentParams,
                                                thread->model);
184
185
186
187

        thread->currentLikelihood = run_state->likelihood(thread->currentParams,
                                                          run_state->data,
                                                          thread->model);
Ben Farr's avatar
Ben Farr committed
188
189
190

        LALInferenceClearVariables(priorDraw);
        XLALFree(priorDraw);
191
192
193
194
    }

    /* Replace distance prior if changed for initial sample draw */
    if (changed_dist) {
Ben Farr's avatar
Ben Farr committed
195
196
        LALInferenceRemoveMinMaxPrior(run_state->priorArgs, "logdistance");
        LALInferenceAddMinMaxPrior(run_state->priorArgs, "logdistance",
197
198
199
                                   &dist_low, &dist_high, LALINFERENCE_REAL8_t);
    }
}
200

201
/*
202
203
204
 * Initialize threads in memory, using LALInferenceInitCBCModel() to init models.
 */
void LALInferenceInitCBCThreads(LALInferenceRunState *run_state, INT4 nthreads) {
Vivien Raymond's avatar
Vivien Raymond committed
205
206
  if (run_state == NULL){
    LALInferenceInitCBCModel(run_state);
Vivien Raymond's avatar
Vivien Raymond committed
207
    return;
Vivien Raymond's avatar
Vivien Raymond committed
208
  }
209
210
211

  ProcessParamsTable *commandLine=run_state->commandLine;

212
  LALInferenceThreadState *thread;
Ben Farr's avatar
Ben Farr committed
213
  INT4 t, nifo;
214
  INT4 randomseed;
215
  LALInferenceIFOData *data;
216
217
  run_state->nthreads = nthreads;
  run_state->threads = LALInferenceInitThreads(nthreads);
218

219
  for (t = 0; t < nthreads; t++) {
220
    thread = &(run_state->threads[t]);
221

222
223
224
    /* Link back to run-state */
    thread->parent = run_state;

225
226
    /* Set up CBC model and parameter array */
    thread->model = LALInferenceInitCBCModel(run_state);
227
    thread->model->roq_flag = 0;
228

Ben Farr's avatar
Ben Farr committed
229
230
    /* Allocate IFO likelihood holders */
    nifo = 0;
231
    data = run_state->data;
Ben Farr's avatar
Ben Farr committed
232
233
234
235
236
237
    while (data != NULL) {
        data = data->next;
        nifo++;
    }
    thread->currentIFOLikelihoods = XLALCalloc(nifo, sizeof(REAL8));

238
    /* Setup ROQ */
239
240
241
242
243
    if (LALInferenceGetProcParamVal(commandLine, "--roqtime_steps")){

        LALInferenceSetupROQmodel(thread->model, commandLine);
        fprintf(stderr, "done LALInferenceSetupROQmodel\n");

244
245
246
    }else{
      thread->model->roq_flag=0;
    }
247

248
    LALInferenceCopyVariables(thread->model->params, thread->currentParams);
249
    LALInferenceCopyVariables(run_state->proposalArgs, thread->proposalArgs);
250

251
252
253
    /* Link thread-state prior-args to the parent run-state's */
    thread->priorArgs = run_state->priorArgs;

254
    /* Use clocktime if seed isn't provided */
255
256
257
    thread->GSLrandom = gsl_rng_alloc(gsl_rng_mt19937);
    randomseed = gsl_rng_get(run_state->GSLrandom);
    gsl_rng_set(thread->GSLrandom, randomseed);
258
  }
259

260
  return;
261
262
263
}


264
265
/* Setup the template generation */
/* Defaults to using LALSimulation */
266
LALInferenceTemplateFunction LALInferenceInitCBCTemplate(LALInferenceRunState *runState)
267
{
268
  char help[]="(--template [LAL,PhenSpin,LALGenerateInspiral,LALSim,multiband]\tSpecify template (default LAL)\n";
269
270
271
272
273
274
275
276
277
278
279
  ProcessParamsTable *ppt=NULL;
  ProcessParamsTable *commandLine=runState->commandLine;
  /* Print command line arguments if help requested */
  //Help is taken care of in LALInferenceInitCBCVariables
  //ppt=LALInferenceGetProcParamVal(commandLine,"--help");
  //if(ppt)
  //{
  //	fprintf(stdout,"%s",help);
  //	return;
  //}
  /* This is the LAL template generator for inspiral signals */
280
  LALInferenceTemplateFunction templt = &LALInferenceTemplateXLALSimInspiralChooseWaveform;
281
282
  ppt=LALInferenceGetProcParamVal(commandLine,"--template");
  if(ppt) {
283
    if(!strcmp("LALSim",ppt->value))
284
      templt=&LALInferenceTemplateXLALSimInspiralChooseWaveform;
285
    else if(!strcmp("null",ppt->value))
286
        templt=&LALInferenceTemplateNullFreqdomain;
287
288
289
290
291
	else if(!strcmp("multiband",ppt->value)){
        templt=&LALInferenceTemplateXLALSimInspiralChooseWaveformPhaseInterpolated;
        fprintf(stdout,"Template function called is \"LALInferenceTemplateXLALSimInspiralChooseWaveformPhaseInterpolated\"\n");
    }
    else {
292
293
294
        XLALPrintError("Error: unknown template %s\n",ppt->value);
        XLALPrintError("%s", help);
        XLAL_ERROR_NULL(XLAL_EINVAL);
295
    }
296
  }
297
298
299
300
301
  else if(LALInferenceGetProcParamVal(commandLine,"--LALSimulation")){
    fprintf(stderr,"Warning: --LALSimulation is deprecated, the LALSimulation package is now the default. To use LALInspiral specify:\n\
                    --template LALGenerateInspiral (for time-domain templates)\n\
                    --template LAL (for frequency-domain templates)\n");
  }
302
  else if(LALInferenceGetProcParamVal(commandLine,"--roqtime_steps")){
303
  templt=&LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence;
304
        fprintf(stderr, "template is \"LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence\"\n");
305
  }
306
307
308
  else {
    fprintf(stdout,"Template function called is \"LALInferenceTemplateXLALSimInspiralChooseWaveform\"\n");
  }
309
  return templt;
310
311
}

312
/* Setup the glitch model */
313
void LALInferenceInitGlitchVariables(LALInferenceRunState *runState, LALInferenceVariables *currentParams)
314
315
316
317
318
319
320
321
322
{
  ProcessParamsTable    *commandLine   = runState->commandLine;
  LALInferenceIFOData   *dataPtr       = runState->data;
  LALInferenceVariables *priorArgs     = runState->priorArgs;

  UINT4 i,nifo;
  UINT4 n = (UINT4)dataPtr->timeData->data->length;
  UINT4 gflag  = 1;
  REAL8 gmin   = 0.0;
323
  REAL8 gmax   = 20.0;
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354

  //over-ride default gmax from command line
  if(LALInferenceGetProcParamVal(commandLine, "--glitchNmax"))
    gmax = (REAL8)atoi(LALInferenceGetProcParamVal(commandLine, "--glitchNmax")->value);

  //count interferometers in network before allocating memory
  //compute imin,imax for each IFO -- may be different
  nifo=0;
  dataPtr = runState->data;
  while (dataPtr != NULL)
  {
    dataPtr = dataPtr->next;
    nifo++;
  }
  dataPtr = runState->data;

  UINT4Vector *gsize  = XLALCreateUINT4Vector(nifo);
  //Meyer?? REAL8Vector *gprior = XLALCreateREAL8Vector((int)gmax+1);

  //Morlet??
  gsl_matrix *mAmp = gsl_matrix_alloc(nifo,(int)(gmax));
  gsl_matrix *mf0  = gsl_matrix_alloc(nifo,(int)(gmax));
  gsl_matrix *mQ   = gsl_matrix_alloc(nifo,(int)(gmax));
  gsl_matrix *mt0  = gsl_matrix_alloc(nifo,(int)(gmax));
  gsl_matrix *mphi = gsl_matrix_alloc(nifo,(int)(gmax));

  double Amin,Amax;
  double Qmin,Qmax;
  double f_min,f_max;
  double tmin,tmax;
  double pmin,pmax;
355
  double Anorm;
356
357

  REAL8 TwoDeltaToverN = 2.0 * dataPtr->timeData->deltaT / ((double) dataPtr->timeData->data->length);
358
359
360
361

  Anorm = sqrt(TwoDeltaToverN);
  Amin = 10.0/Anorm;
  Amax = 10000.0/Anorm;
362
363
364
365
366
367
368
369
370

  Qmin = 3.0;
  Qmax = 30.0;
  tmin = 0.0;
  tmax = dataPtr->timeData->data->length*dataPtr->timeData->deltaT;
  f_min = dataPtr->fLow;
  f_max = dataPtr->fHigh;
  pmin = 0.0;
  pmax = LAL_TWOPI;
371

372
373
374
375
376
377
  gsl_matrix_set_all(mAmp, Amin);
  gsl_matrix_set_all(mf0,  f_min);
  gsl_matrix_set_all(mQ,   Qmin);
  gsl_matrix_set_all(mt0,  tmin);
  gsl_matrix_set_all(mphi, pmin);

378
379
380
381
382
383
384
385
386
387
388
389
  gsl_matrix  *gFD       = gsl_matrix_alloc(nifo,(int)n); //store the Fourier-domain glitch signal

  for(i=0; i<nifo; i++) gsize->data[i]=0;

  //Morlet wavelet parameters
  LALInferenceAddVariable(currentParams, "morlet_FD",  &gFD,  LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);
  LALInferenceAddVariable(currentParams, "morlet_Amp", &mAmp, LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);
  LALInferenceAddVariable(currentParams, "morlet_f0" , &mf0,  LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);
  LALInferenceAddVariable(currentParams, "morlet_Q"  , &mQ,   LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);
  LALInferenceAddVariable(currentParams, "morlet_t0" , &mt0,  LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);
  LALInferenceAddVariable(currentParams, "morlet_phi", &mphi, LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);

390
  LALInferenceAddVariable(currentParams, "glitch_size", &gsize, LALINFERENCE_UINT4Vector_t, LALINFERENCE_PARAM_LINEAR);
391
392
393
394
395
396
397
398
  LALInferenceAddVariable(currentParams, "glitchFitFlag", &gflag, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);

  LALInferenceAddMinMaxPrior(priorArgs, "morlet_Amp_prior", &Amin, &Amax, LALINFERENCE_REAL8_t);
  LALInferenceAddMinMaxPrior(priorArgs, "morlet_f0_prior" , &f_min, &f_max, LALINFERENCE_REAL8_t);
  LALInferenceAddMinMaxPrior(priorArgs, "morlet_Q_prior"  , &Qmin, &Qmax, LALINFERENCE_REAL8_t);
  LALInferenceAddMinMaxPrior(priorArgs, "morlet_t0_prior" , &tmin, &tmax, LALINFERENCE_REAL8_t);
  LALInferenceAddMinMaxPrior(priorArgs, "morlet_phi_prior", &pmin, &pmax, LALINFERENCE_REAL8_t);

399
  LALInferenceAddMinMaxPrior(priorArgs, "glitch_size", &gmin, &gmax, LALINFERENCE_REAL8_t);
400
  LALInferenceAddMinMaxPrior(priorArgs, "glitch_dim", &gmin, &gmax, LALINFERENCE_REAL8_t);
401

402
  LALInferenceAddVariable(priorArgs, "glitch_norm", &Anorm, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
403
404
}

405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
struct spcal_envelope
{
    gsl_spline  *amp_median,*amp_std,
                *phase_median,*phase_std;
};

/* Format string for the calibratino envelope file */
/* Frequency    Median Mag     Phase (Rad)    -1 Sigma Mag   -1 Sigma Phase +1 Sigma Mag   +1 Sigma Phase */

#define CAL_ENV_FORMAT "%lf %lf %lf %lf %lf %lf %lf\n"

static struct spcal_envelope *initCalibrationEnvelope(char *filename);

static struct spcal_envelope *initCalibrationEnvelope(char *filename)
{
    FILE *fp=fopen(filename,"r");
    char tmpline[1024];
    if(!fp) {fprintf(stderr,"Unable to open %s: Error %i %s\n",filename,errno,strerror(errno)); exit(1);}
    int Nlines=0;
    REAL8 freq, *logfreq=NULL, *mag_med=NULL, mag_low, mag_hi, *mag_std=NULL, *phase_med=NULL, phase_low, phase_hi, *phase_std=NULL;
    for(Nlines=0;fgets(tmpline,1024,fp); )
    {
        /* Skip header */
        if(tmpline[0]=='#') continue;
        /* Grow arrays */
        logfreq=realloc(logfreq,sizeof(*logfreq)*(Nlines+1));
        mag_med=realloc(mag_med,sizeof(*mag_med)*(Nlines+1));
        mag_std=realloc(mag_std,sizeof(*mag_std)*(Nlines+1));
        phase_med=realloc(phase_med,sizeof(*phase_med)*(Nlines+1));
        phase_std=realloc(phase_std,sizeof(*phase_std)*(Nlines+1));

        if((7!=sscanf(tmpline,CAL_ENV_FORMAT, &freq, &(mag_med[Nlines]), &(phase_med[Nlines]), &mag_low, &phase_low, &mag_hi, &phase_hi)))
        {
            fprintf(stderr,"Malformed input line in file %s: %s\n",filename,tmpline);
            exit(1);
        }
441
		mag_med[Nlines]-=1.0; /* Subtract off 1 to get delta */
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
        logfreq[Nlines]=log(freq);
        mag_std[Nlines]=(mag_hi - mag_low ) /2.0;
        phase_std[Nlines]=(phase_hi - phase_low) /2.0;
		Nlines++;
    }
    fprintf(stdout,"Read %i lines from calibration envelope %s\n",Nlines,filename);
    fclose(fp);

    struct spcal_envelope *env=XLALCalloc(1,sizeof(*env));
    env->amp_median = gsl_spline_alloc ( gsl_interp_cspline, Nlines);
    env->amp_std = gsl_spline_alloc ( gsl_interp_cspline, Nlines);
    env->phase_median = gsl_spline_alloc ( gsl_interp_cspline, Nlines);
    env->phase_std = gsl_spline_alloc ( gsl_interp_cspline, Nlines);

    gsl_spline_init(env->amp_median, logfreq, mag_med, Nlines);
    gsl_spline_init(env->amp_std, logfreq, mag_std, Nlines);
    gsl_spline_init(env->phase_median, logfreq, phase_med, Nlines);
    gsl_spline_init(env->phase_std, logfreq, phase_std, Nlines);
460

461
    free(logfreq); free(mag_med); free(mag_std); free(phase_med); free(phase_std);
462

463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
    return(env);
}

static int destroyCalibrationEnvelope(struct spcal_envelope *env);
static int destroyCalibrationEnvelope(struct spcal_envelope *env)
{
    if(!env) XLAL_ERROR(XLAL_EINVAL);
    if(env->amp_median) gsl_spline_free(env->amp_median);
    if(env->amp_std) gsl_spline_free(env->amp_std);
    if(env->phase_median) gsl_spline_free(env->phase_median);
    if(env->phase_std) gsl_spline_free(env->phase_std);
    XLALFree(env);
    return XLAL_SUCCESS;
}

478
void LALInferenceInitCalibrationVariables(LALInferenceRunState *runState, LALInferenceVariables *currentParams) {
479
480
  ProcessParamsTable *ppt = NULL;
  LALInferenceIFOData *ifo = NULL;
481
482
483
484
485
486
487
488
  LALInferenceIFOData *dataPtr=NULL;
  UINT4 calOn = 1;
  if ((ppt = LALInferenceGetProcParamVal(runState->commandLine, "--enable-spline-calibration"))){
    /* Use spline to marginalize*/
    UINT4 ncal = 5; /* Number of calibration nodes, log-distributed
		between fmin and fmax. */
    REAL8 ampUncertaintyPrior = 0.1; /* 10% amplitude */
    REAL8 phaseUncertaintyPrior = 5*M_PI/180.0; /* 5 degrees phase */
489
    if ((ppt = LALInferenceGetProcParamVal(runState->commandLine, "--spcal-nodes"))) {
490
      ncal = atoi(ppt->value);
491
492
493
494
      if (ncal < 3) { /* Cannot do spline with fewer than 3 points! */
	fprintf(stderr, "ERROR: given '--spcal-nodes %d', but cannot spline with fewer than 3\n", ncal);
	exit(1);
      }
495
    }
496

497
498
499
    LALInferenceAddVariable(currentParams, "spcal_active", &calOn, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);
    LALInferenceAddVariable(currentParams, "spcal_npts", &ncal, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);

500
    for(ifo=runState->data;ifo;ifo=ifo->next) {
501
      UINT4 i;
502
503
504
505
506
507
508
509
510
511
512

      char freqVarName[VARNAME_MAX];
      char ampVarName[VARNAME_MAX];
      char phaseVarName[VARNAME_MAX];

      REAL8 fMin = ifo->fLow;
      REAL8 fMax = ifo->fHigh;
      REAL8 logFMin = log(fMin);
      REAL8 logFMax = log(fMax);
      REAL8 dLogF = (logFMax - logFMin)/(ncal-1);

513
514
      char amp_uncert_op[VARNAME_MAX];
      char pha_uncert_op[VARNAME_MAX];
515
516
      char env_uncert_op[VARNAME_MAX];
      struct spcal_envelope *env=NULL;
517

518
519
520
521
522
523
      if((VARNAME_MAX <= snprintf(amp_uncert_op, VARNAME_MAX, "--%s-spcal-amp-uncertainty", ifo->name))
          || (VARNAME_MAX <= snprintf(pha_uncert_op, VARNAME_MAX, "--%s-spcal-phase-uncertainty", ifo->name))
          || (VARNAME_MAX <= snprintf(env_uncert_op, VARNAME_MAX, "--%s-spcal-envelope",ifo->name)) )
      {
        fprintf(stderr,"variable name too long\n"); exit(1);
      }
524

525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
      if( (ppt=LALInferenceGetProcParamVal(runState->commandLine, env_uncert_op)))
          env = initCalibrationEnvelope(ppt->value);
      else
      {
        if ((ppt = LALInferenceGetProcParamVal(runState->commandLine, amp_uncert_op))) {
            ampUncertaintyPrior = atof(ppt->value);
        }
        else{
            fprintf(stderr,"Error, missing %s or %s\n",amp_uncert_op, env_uncert_op);
            exit(1);
        }

        if ((ppt = LALInferenceGetProcParamVal(runState->commandLine, pha_uncert_op))) {
            phaseUncertaintyPrior = M_PI/180.0*atof(ppt->value); /* CL arg in degrees, variable in radians */
        }
        else{
            fprintf(stderr,"Error, missing %s or %s\n",pha_uncert_op,env_uncert_op);
            exit(1);
        }
544
      }
545
546
      /* Now add each spline node */
      for(i=0;i<ncal;i++)
547
	  {
548
549
550
551
552
553
554
			  if((VARNAME_MAX <= snprintf(freqVarName, VARNAME_MAX, "%s_spcal_logfreq_%i",ifo->name,i))
                  || (VARNAME_MAX <= snprintf(ampVarName, VARNAME_MAX, "%s_spcal_amp_%i", ifo->name,i))
                  || (VARNAME_MAX <= snprintf(phaseVarName, VARNAME_MAX, "%s_spcal_phase_%i", ifo->name,i)) )
              {
                  fprintf(stderr,"Variable name too long\n"); exit(1);
              }
              REAL8 amp_std=ampUncertaintyPrior,amp_mean=0.0;
555
556
557
558
559
560
561
562
			  REAL8 phase_std=phaseUncertaintyPrior,phase_mean=0.0;
			  REAL8 logFreq = logFMin + i*dLogF;
			  LALInferenceAddREAL8Variable(currentParams,freqVarName,logFreq,LALINFERENCE_PARAM_FIXED);
			  if(env)
			  {
					  amp_std = gsl_spline_eval(env->amp_std, logFreq, NULL);
					  amp_mean = gsl_spline_eval(env->amp_median, logFreq, NULL);
					  phase_std = gsl_spline_eval(env->phase_std, logFreq, NULL);
563
					  phase_mean = gsl_spline_eval(env->phase_median, logFreq, NULL);
564
565
566
567
568
569
570
			  }
			  LALInferenceRegisterGaussianVariableREAL8(runState, currentParams, ampVarName, 0, amp_mean, amp_std, LALINFERENCE_PARAM_LINEAR);
			  LALInferenceRegisterGaussianVariableREAL8(runState, currentParams, phaseVarName, 0, phase_mean, phase_std, LALINFERENCE_PARAM_LINEAR);
	  } /* End loop over spline nodes */

	  if(env) destroyCalibrationEnvelope(env);
	} /* End loop over IFOs */
571
  } /* End case of spline calibration error */
572
573
574
575
576
577
578
579
580
  else if(LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalAmp") ||LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalPha")){
    /* Use constant (in frequency) approximation for the errors */
    if (LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalAmp")){
      /*For the moment the prior ranges are the same for the three IFOs */
      REAL8 camp_max_A=0.25; /* plus minus 25% amplitude errors*/
      REAL8 camp_min_A=-0.25;
      REAL8 zero=0.0;
      dataPtr = runState->data;
      while (dataPtr != NULL){
581
        char CA_A[320];
582
583
584
585
        sprintf(CA_A,"%s_%s","calamp",dataPtr->name);
        LALInferenceRegisterUniformVariableREAL8(runState, currentParams, CA_A, zero, camp_min_A, camp_max_A, LALINFERENCE_PARAM_LINEAR);
        dataPtr = dataPtr->next;
      }
586

587
      LALInferenceAddVariable(currentParams, "constantcal_active", &calOn, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);
588
589
590
591
592
593
594
      /*If user specifies a width for the error prior, a gaussian prior will be used, otherwise a flat prior will be used*/
      REAL8 ampUncertaintyPrior=-1.0;
      ppt = LALInferenceGetProcParamVal(runState->commandLine, "--constcal_ampsigma");
      if (ppt) {
        ampUncertaintyPrior = atof(ppt->value);
      }
      LALInferenceAddVariable(runState->priorArgs, "constcal_amp_uncertainty", &ampUncertaintyPrior, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
595
    }
596
597
598
599
600
601
602
603
    if (LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalPha")){
      /* Add linear calibration phase errors to the measurement. For the moment the prior ranges are the same for the three IFOs */
      REAL8 cpha_max_A=0.349;  /* plus/minus 20 degs*/
      REAL8 cpha_min_A=-0.349;
      REAL8 zero=0.0;
      dataPtr = runState->data;
      while (dataPtr != NULL)
      {
604
        char CP_A[320];
605
606
607
608
609
        sprintf(CP_A,"%s_%s","calpha",dataPtr->name);
        LALInferenceRegisterUniformVariableREAL8(runState, currentParams, CP_A, zero, cpha_min_A, cpha_max_A, LALINFERENCE_PARAM_LINEAR);
        dataPtr = dataPtr->next;
      }
      LALInferenceAddVariable(currentParams, "constantcal_active", &calOn, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);
610
611
612
613
614
615
616
617
618

     /*If user specifies a width for the error prior, a gaussian prior will be used, otherwise a flat prior will be used*/
      REAL8 phaseUncertaintyPrior=-1.0;
      ppt = LALInferenceGetProcParamVal(runState->commandLine, "--constcal_phasigma");
      if (ppt) {
        phaseUncertaintyPrior = M_PI/180.0*atof(ppt->value); /* CL arg in degrees, variable in radians */
      }
      LALInferenceAddVariable(runState->priorArgs, "constcal_phase_uncertainty", &phaseUncertaintyPrior, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);

619
    }
620
621
622
623
624
  }
  else{
    /* No calibration marginalization asked. Just exit */
    return;
  }
625
626
}

627
628
629
630
631
632
633
void LALInferenceRegisterGaussianVariableREAL8(LALInferenceRunState *state, LALInferenceVariables *var, const char name[VARNAME_MAX], REAL8 startval, REAL8 mean, REAL8 stdev, LALInferenceParamVaryType varytype)
{
  char meanopt[VARNAME_MAX+8];
  char sigmaopt[VARNAME_MAX+9];
  char valopt[VARNAME_MAX+3];
  char fixopt[VARNAME_MAX+7];
  ProcessParamsTable *ppt=NULL;
634

635
636
637
638
  sprintf(meanopt,"--%s-mean",name);
  sprintf(sigmaopt,"--%s-sigma",name);
  sprintf(valopt,"--%s",name);
  sprintf(fixopt,"--fix-%s",name);
639

640
641
642
643
644
645
646
647
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,meanopt))) mean=atof(ppt->value);
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,sigmaopt))) stdev=atof(ppt->value);
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,fixopt)))
  {
    varytype = LALINFERENCE_PARAM_FIXED;
    startval = atof(ppt->value);
  }
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,valopt))) startval=atof(ppt->value);
648

649
650
651
652
653
654
  assert(stdev>0);
  LALInferenceAddVariable(var,name,&startval,LALINFERENCE_REAL8_t,varytype);
  LALInferenceAddGaussianPrior(state->priorArgs, name, &mean, &stdev, LALINFERENCE_REAL8_t);

}

655
656
657
658
659
660
661
void LALInferenceRegisterUniformVariableREAL8(LALInferenceRunState *state, LALInferenceVariables *var, const char name[VARNAME_MAX], REAL8 startval, REAL8 min, REAL8 max, LALInferenceParamVaryType varytype)
{
  char minopt[VARNAME_MAX+7];
  char maxopt[VARNAME_MAX+7];
  char valopt[VARNAME_MAX+3];
  char fixopt[VARNAME_MAX+7];
  ProcessParamsTable *ppt=NULL;
662

663
664
665
666
  sprintf(minopt,"--%s-min",name);
  sprintf(maxopt,"--%s-max",name);
  sprintf(valopt,"--%s",name);
  sprintf(fixopt,"--fix-%s",name);
667

668
669
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,minopt))) min=atof(ppt->value);
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,maxopt))) max=atof(ppt->value);
670
671
672
673
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,fixopt))) {
		  varytype=LALINFERENCE_PARAM_FIXED;
		  startval=atof(ppt->value);
  }
674
675
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,valopt))) startval=atof(ppt->value);
  else if(varytype!=LALINFERENCE_PARAM_FIXED) startval=min+(max-min)*gsl_rng_uniform(state->GSLrandom);
676

677
678
679
680
681
682
683
684
685
  /* Error checking */
  if(min>max) {
    fprintf(stderr,"ERROR: Prior for %s has min(%lf) > max(%lf)\n",name,min,max);
    exit(1);
  }
  if(startval<min || startval>max){
    fprintf(stderr,"ERROR: Initial value %lf for %s lies outwith prior (%lf,%lf)\n",startval,name,min,max);
    exit(1);
  }
686
687
688
689
690
691
692
693
694
  /* Mass parameters checks*/
  if (!strcmp(name,"eta"))
    if (max>0.25){
      fprintf(stderr,"ERROR: maximum of eta cannot be larger than 0.25. Check --eta-max\n");
      exit(1);
    }
  if (!strcmp(name,"q")){
    REAL8 qMin=min;
    REAL8 qMax=max;
695

696
    if (qMin <= 0.0 || qMin > 1.0)
697
    {
698
699
        fprintf(stderr,"ERROR: qMin must be between 0 and 1, got value qMin=%f\n",qMin);
		exit(1);
700
    }
701
    if (qMax > 1.0 || qMax <0.0 || qMax < qMin)
702
    {
703
704
      fprintf(stderr,"ERROR: qMax must be between 0 and 1, and qMax > qMin. Got value qMax=%f, qMin=%f\n",qMax,qMin);
	  exit(1);
705
706
707
708
    }
  }
  /*End of mass parameters check */

709
710
  LALInferenceAddVariable(var,name,&startval,LALINFERENCE_REAL8_t,varytype);
  LALInferenceAddMinMaxPrior(state->priorArgs, name, &min, &max, LALINFERENCE_REAL8_t);
711

712
713
}

714

715
/* Setup the variables to control template generation for the CBC model */
716
/* Includes specification of prior ranges. Returns address of new LALInferenceVariables */
717

Ben Farr's avatar
Ben Farr committed
718
LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
719
  char help[]="\
Vivien Raymond's avatar
Vivien Raymond committed
720
721
722
723
724
    ----------------------------------------------\n\
    --- Injection Arguments ----------------------\n\
    ----------------------------------------------\n\
    (--inj injections.xml) Injection XML file to use\n\
    (--event N)            Event number from Injection XML file to use\n\
725
\n\
Vivien Raymond's avatar
Vivien Raymond committed
726
727
728
729
730
731
    ----------------------------------------------\n\
    --- Template Arguments -----------------------\n\
    ----------------------------------------------\n\
    (--use-eta)            Jump in symmetric mass ratio eta, instead of q=m1/m2 (m1>m2)\n\
    (--approx)             Specify a template approximant and phase order to use\n\
                         (default TaylorF2threePointFivePN). Available approximants:\n\
732
733
734
735
736
737
738
739
                         default modeldomain=\"time\": GeneratePPN, TaylorT1, TaylorT2,\n\
                                                       TaylorT3, TaylorT4, EOB, EOBNR,\n\
                                                       EOBNRv2, EOBNRv2HM, SEOBNRv1,\n\
                                                       SpinTaylor, SpinQuadTaylor, \n\
                                                       SpinTaylorFrameless, SpinTaylorT4,\n\
                                                       PhenSpinTaylorRD, NumRel.\n\
                         default modeldomain=\"frequency\": TaylorF1, TaylorF2, TaylorF2RedSpin,\n\
                                                       TaylorF2RedSpinTidal, IMRPhenomA,\n\
740
                                                       IMRPhenomB, IMRPhenomP, IMRPhenomPv2.\n\
Vivien Raymond's avatar
Vivien Raymond committed
741
742
    (--amporder PNorder)            Specify a PN order in amplitude to use (defaults: LALSimulation: max available; LALInspiral: newtownian).\n\
    (--fref f_ref)                  Specify a reference frequency at which parameters are defined (default 100).\n\
Salvatore Vitale's avatar
Salvatore Vitale committed
743
744
    (--tidal)                   Enables tidal corrections, only with LALSimulation.\n\
    (--tidalT)                  Enables reparmeterized tidal corrections, only with LALSimulation.\n\
745
    (--4PolyEOS)                Enables 4-piece polytropic EOS parmeterization, only with LALSimulation.\n\
John Douglas Veitch's avatar
John Douglas Veitch committed
746
747
    (--4SpectralDecomp)         Enables 4-coeff. spectral decomposition EOS parmeterization, only with LALSimulation.\n\
    (--eos   EOS)               Fix the neutron star EOS. Use \"--eos help\" for allowed names\n\
Carl-Johan Haster's avatar
Carl-Johan Haster committed
748
    (--BinaryLove)              Enable the Binary Neutron Star common EoS tidal model, expressed through EOS-insensitive relations\n\
Vivien Raymond's avatar
Vivien Raymond committed
749
750
    (--spinOrder PNorder)           Specify twice the PN order (e.g. 5 <==> 2.5PN) of spin effects to use, only for LALSimulation (default: -1 <==> Use all spin effects).\n\
    (--tidalOrder PNorder)          Specify twice the PN order (e.g. 10 <==> 5PN) of tidal effects to use, only for LALSimulation (default: -1 <==> Use all tidal effects).\n\
751
    (--numreldata FileName)         Location of NR data file for NR waveforms (with NR_hdf5 approx).\n\
Vivien Raymond's avatar
Vivien Raymond committed
752
753
754
755
756
    (--modeldomain)                 domain the waveform template will be computed in (\"time\" or \"frequency\"). If not given will use LALSim to decide\n\
    (--spinAligned or --aligned-spin)  template will assume spins aligned with the orbital angular momentum.\n\
    (--singleSpin)                  template will assume only the spin of the most massive binary component exists.\n\
    (--noSpin, --disable-spin)      template will assume no spins (giving this will void spinOrder!=0) \n\
    (--no-detector-frame)              model will NOT use detector-centred coordinates and instead RA,dec\n\
757
758
759
    (--nonGR_alpha value) this is a LIV parameter which should only be passed when log10lambda_eff/lambda_eff is passed as a grtest-parameter for LIV test\n\
    (--LIV_A_sign) this is a LIV parameter determining if +A or -A is being tested; A occurs in the modified dispersion relation. LIV_A_sign has to be either +1 or -1 \n\
    (--liv) this flag is now needed for launching a LIV run\n\
John Douglas Veitch's avatar
John Douglas Veitch committed
760
761
    (--grtest-parameters dchi0,..,dxi1,..,dalpha1,..) template will assume deformations in the corresponding phase coefficients.\n\
    (--ppe-parameters aPPE1,....     template will assume the presence of an arbitrary number of PPE parameters. They must be paired correctly.\n\
762
    (--modeList lm,l-m...,lm,l-m)           List of modes to be used by the model. The chosen modes ('lm') should be passed as a ',' seperated list.\n\
763
764
765
766
767
768
769
770
    (--phenomXHMMband float)     Threshold parameter for the Multibanding of the non-precessing hlm modes in IMRPhenomXHM and IMRPhenomXPHM. If set to 0 then do not use multibanding.\n\
                                 Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
    (--phenomXPHMMband float)    Threshold parameter for the Multibanding of the Euler angles in IMRPhenomXPHM. If set to 0 then do not use multibanding.\n\
                                 Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
    (--phenomXPFinalSpinMod int) Change version for the final spin model used in IMRPhenomXP/IMRPhenomXPHM.\n\
                                 Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
    (--phenomXPrecVersion int)   Change version of the Euler angles for the twisting-up of IMRPhenomXP/IMRPhenomXPHM.\n\
                                 Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
771
\n\
Vivien Raymond's avatar
Vivien Raymond committed
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
    ----------------------------------------------\n\
    --- Starting Parameters ----------------------\n\
    ----------------------------------------------\n\
    You can generally have MCMC chains to start from a given parameter value by using --parname VALUE. Names currently known to the code are:\n\
     time                         Waveform time (overrides random about trigtime).\n\
     chirpmass                    Chirpmass\n\
     eta                          Symmetric massratio (needs --use-eta)\n\
     q                            Asymmetric massratio (a.k.a. q=m2/m1 with m1>m2)\n\
     phase                        Coalescence phase.\n\
     costheta_jn                  Cosine of angle between J and line of sight [rads]\n\
     logdistance                  Log Distance (requires --use-logdistance)\n\
     rightascension               Rightascensions\n\
     declination                  Declination.\n\
     polarisation                 Polarisation angle.\n\
    * Spin Parameters:\n\
     a_spin1                      Spin1 magnitude\n\
     a_spin2                      Spin2 magnitude\n\
     tilt_spin1                   Angle between spin1 and orbital angular momentum\n\
     tilt_spin2                   Angle between spin2 and orbital angular momentum \n\
     phi_12                       Difference between spins' azimuthal angles \n\
     phi_jl                       Difference between total and orbital angular momentum azimuthal angles\n\
793
    * Equation of State parameters:\n\
794
     (requires --tidal)\n\
Vivien Raymond's avatar
Vivien Raymond committed
795
796
     lambda1                      lambda1.\n\
     lambda2                      lambda2.\n\
797
     (requires --tidalT)\n\
Vivien Raymond's avatar
Vivien Raymond committed
798
799
     lambdaT                      lambdaT.\n\
     dLambdaT                     dLambdaT.\n\
800
     (requires --4PolyEOS)\n\
801
802
803
804
     logp1                        logp1.\n\
     gamma1                       gamma1.\n\
     gamma2                       gamma2.\n\
     gamma3                       gamma3.\n\
Carl-Johan Haster's avatar
Carl-Johan Haster committed
805
806
     (requires --BinaryLove)\n\
     lambdaS                      Symmetric tidal deformability.\n\
807
808
809
810
811
     (requires --4SpectralDecomp):\n\
     SDgamma0                     SDgamma0.\n\
     SDgamma1                     SDgamma1.\n\
     SDgamma2                     SDgamma2.\n\
     SDgamma3                     SDgamma3.\n\
812
    * \n\
NV KRISHNENDU's avatar
NV KRISHNENDU committed
813
814
815
816
817
818
819
820
    * Spin-induced quadrupole moment test parameters:\n\
     (requires --dQuadMon12)\n\
     dQuadMon1                      dQuadMon1.\n\
     dQuadMon2                      dQuadMon2.\n\
     (requires --dQuadMonSA) \n\
     dQuadMonS                      dQuadMonS.\n\
     dQuadMonA                      dQuadMonA.\n\
    (dQuadMonS and dQuadMonA are the symmetric and antisymmetric combinations of dQuadMon1 and dQuadMon2).\n\
Vivien Raymond's avatar
Vivien Raymond committed
821
822
823
824
825
826
827
828
    ----------------------------------------------\n\
    --- Prior Ranges -----------------------------\n\
    ----------------------------------------------\n\
    You can generally use --paramname-min MIN --paramname-max MAX to set the prior range for the parameter paramname\n\
    The names known to the code are listed below.\n\
    Component masses, total mass and time have dedicated options listed here:\n\n\
    (--trigtime time)                       Center of the prior for the time variable.\n\
    (--comp-min min)                        Minimum component mass (1.0).\n\
829
    (--comp-max max)                        Maximum component mass (100.0).\n\
Vivien Raymond's avatar
Vivien Raymond committed
830
831
832
    (--mass1-min min, --mass1-max max)      Min and max for mass1 (default: same as comp-min,comp-max, will over-ride these.\n\
    (--mass2-min min, --mass2-max max)      Min and max for mass2 (default: same as comp-min,comp-max, will over-ride these.\n\
    (--mtotal-min min)                      Minimum total mass (2.0).\n\
833
    (--mtotal-max max)                      Maximum total mass (200.0).\n\
Vivien Raymond's avatar
Vivien Raymond committed
834
    (--dt time)                             Width of time prior, centred around trigger (0.2s).\n\
835
\n\
836
    (--ns-max-mass)                          Maximum observed NS mass used for EOS prior (none).\n\
Vivien Raymond's avatar
Vivien Raymond committed
837
838
839
840
841
842
843
    (--varyFlow, --flowMin, --flowMax)       Allow the lower frequency bound of integration to vary in given range.\n\
    (--pinparams)                            List of parameters to set to injected values [mchirp,asym_massratio,etc].\n\
    ----------------------------------------------\n\
    --- Fix Parameters ---------------------------\n\
    ----------------------------------------------\n\
    You can generally fix a parameter to be fixed to a given values by using --fix-paramname VALUE\n\
    where the known names have been listed above.\n\
844
\n";
845
846
847
848
849

  /* Print command line arguments if state was not allocated */
  if(state==NULL)
    {
      fprintf(stdout,"%s",help);
850
      return(NULL);
851
852
853
854
855
856
    }

  /* Print command line arguments if help requested */
  if(LALInferenceGetProcParamVal(state->commandLine,"--help"))
    {
      fprintf(stdout,"%s",help);
857
      return(NULL);
858
859
860
861
    }

  LALStatus status;
  memset(&status,0,sizeof(status));
862
  int errnum;
863
864
  SimInspiralTable *injTable=NULL;
  LALInferenceVariables *priorArgs=state->priorArgs;
865
  LALInferenceVariables *proposalArgs=state->proposalArgs;
866
867
868
  ProcessParamsTable *commandLine=state->commandLine;
  ProcessParamsTable *ppt=NULL;
  ProcessParamsTable *ppt_order=NULL;
869
870
  LALPNOrder PhaseOrder=-1;
  LALPNOrder AmpOrder=-1;
871
  Approximant approx=NumApproximants;
872
  REAL8 f_ref = 100.0;
873
874
875
  LALInferenceIFOData *dataPtr;
  UINT4 event=0;
  UINT4 i=0;
876
  /* Default priors */
877
  REAL8 Dmin=1.0;
878
  REAL8 Dmax=2000.0;
879
  REAL8 Dinitial = (Dmax + Dmin)/2.0;
880
881
882
883
  REAL8 mcMin=1.0;
  REAL8 mcMax=15.3;
  REAL8 etaMin=0.0312;
  REAL8 etaMax=0.25;
884
  REAL8 qMin=1./30.; // The ratio between min and max component mass (see InitMassVariables)
885
886
  REAL8 qMax=1.0;
  REAL8 psiMin=0.0,psiMax=LAL_PI;
John Douglas Veitch's avatar
John Douglas Veitch committed
887
888
  REAL8 decMin=-LAL_PI/2.0,decMax=LAL_PI/2.0;
  REAL8 raMin=0.0,raMax=LAL_TWOPI;
889
  REAL8 phiMin=0.0,phiMax=LAL_TWOPI;
890
  REAL8 costhetaJNmin=-1.0 , costhetaJNmax=1.0;
891
  REAL8 dt=0.1;  /* Half the width of time prior */
892
893
894
  REAL8 lambda1Min=0.0;
  REAL8 lambda1Max=3000.0;
  REAL8 lambda2Min=0.0;
895
896
897
898
899
  REAL8 lambda2Max=3000.0;
  REAL8 lambdaTMin=0.0;
  REAL8 lambdaTMax=3000.0;
  REAL8 dLambdaTMin=-500.0;
  REAL8 dLambdaTMax=500.0;
NV KRISHNENDU's avatar
NV KRISHNENDU committed
900
901
  REAL8 dQuadMonMin=-200.0;
  REAL8 dQuadMonMax=200.0;
902
903
904
905
906
907
908
909
  REAL8 logp1Min=33.6;
  REAL8 logp1Max=35.4;
  REAL8 gamma1Min=2.0;
  REAL8 gamma1Max=4.5;
  REAL8 gamma2Min=1.1;
  REAL8 gamma2Max=4.5;
  REAL8 gamma3Min=1.1;
  REAL8 gamma3Max=4.5;
Carl-Johan Haster's avatar
Carl-Johan Haster committed
910
911
  REAL8 lambdaSMin=0.0;
  REAL8 lambdaSMax=5000.0;
912
913
914
915
916
917
918
919
  REAL8 SDgamma0Min=0.2;
  REAL8 SDgamma0Max=2.0;
  REAL8 SDgamma1Min=-1.6;
  REAL8 SDgamma1Max=1.7;
  REAL8 SDgamma2Min=-0.6;
  REAL8 SDgamma2Max=0.6;
  REAL8 SDgamma3Min=-0.02;
  REAL8 SDgamma3Max=0.02;
920
  gsl_rng *GSLrandom=state->GSLrandom;
921
922
  REAL8 endtime=0.0, timeParam=0.0;
  REAL8 timeMin=endtime-dt,timeMax=endtime+dt;
923
  REAL8 zero=0.0; /* just a number that will be overwritten anyway*/
924

925
  /* Over-ride prior bounds if analytic test */
926
927
  if (LALInferenceGetProcParamVal(commandLine, "--correlatedGaussianLikelihood"))
  {
928
    return(LALInferenceInitModelReviewEvidence(state));
929
930
931
  }
  else if (LALInferenceGetProcParamVal(commandLine, "--bimodalGaussianLikelihood"))
  {
932
    return(LALInferenceInitModelReviewEvidence_bimod(state));
933
934
935
  }
  else if (LALInferenceGetProcParamVal(commandLine, "--rosenbrockLikelihood"))
  {
936
    return(LALInferenceInitModelReviewEvidence(state)); /* CHECKME: Use the default prior for unimodal */
937
938
  }

Ben Farr's avatar
Ben Farr committed
939
940
941
  LALInferenceModel *model = XLALMalloc(sizeof(LALInferenceModel));
  model->params = XLALCalloc(1, sizeof(LALInferenceVariables));
  memset(model->params, 0, sizeof(LALInferenceVariables));
942
  model->eos_fam = NULL;
Ben Farr's avatar
Ben Farr committed
943

944
945
946
947
  UINT4 signal_flag=1;
  ppt = LALInferenceGetProcParamVal(commandLine, "--noiseonly");
  if(ppt)signal_flag=0;
  LALInferenceAddVariable(model->params, "signalModelFlag", &signal_flag,  LALINFERENCE_INT4_t,  LALINFERENCE_PARAM_FIXED);
948

949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
  /* Read injection XML file for parameters if specified */
  ppt=LALInferenceGetProcParamVal(commandLine,"--inj");
  if(ppt){
    SimInspiralTableFromLIGOLw(&injTable,ppt->value,0,0);
    if(!injTable){
      fprintf(stderr,"Unable to open injection file %s\n",ppt->value);
      exit(1);
    }
    ppt=LALInferenceGetProcParamVal(commandLine,"--event");
    if(ppt){
      event= atoi(ppt->value);
      fprintf(stderr,"Reading event %d from file\n",event);
      i=0;
      while(i<event) {i++; injTable=injTable->next;} /* select event */
    }
  }

966
967
  /* See if there are any parameters pinned to injection values */
  if((ppt=LALInferenceGetProcParamVal(commandLine,"--pinparams"))){
968
    char *pinned_params=ppt->value;
969
970
971
972
973
974
975
976
977
978
    LALInferenceVariables tempParams;
    memset(&tempParams,0,sizeof(tempParams));
    char **strings=NULL;
    UINT4 N;
    LALInferenceParseCharacterOptionString(pinned_params,&strings,&N);
    LALInferenceInjectionToVariables(injTable,&tempParams);
    LALInferenceVariableItem *node=NULL;
    while(N>0){
      N--;
      char *name=strings[N];
979
      fprintf(stdout,"Pinning parameter %s\n",node->name);
980
      node=LALInferenceGetItem(&tempParams,name);
981
      if(node) LALInferenceAddVariable(model->params,node->name,node->value,node->type,node->vary);
982
983
984
      else {fprintf(stderr,"Error: Cannot pin parameter %s. No such parameter found in injection!\n",node->name);}
    }
  }
985

986
  /* Over-ride approximant if user specifies */
987
  ppt=LALInferenceGetProcParamVal(commandLine,"--approximant");
988
989
990
991
992
  if(ppt){
    approx = XLALGetApproximantFromString(ppt->value);
    ppt_order=LALInferenceGetProcParamVal(commandLine,"--order");
    if(ppt_order) PhaseOrder = XLALGetOrderFromString(ppt_order->value);
  }
993
994
995
  ppt=LALInferenceGetProcParamVal(commandLine,"--approx");
  if(ppt){
    approx = XLALGetApproximantFromString(ppt->value);
996
    XLAL_TRY(PhaseOrder = XLALGetOrderFromString(ppt->value),errnum);
997
    if( (int) PhaseOrder == XLAL_FAILURE || errnum) {
998
999
      PhaseOrder=-1;
    }
1000
  }
For faster browsing, not all history is shown. View entire blame