LALInference  4.1.6.1-89842e6
LALInferencePrior.c
Go to the documentation of this file.
1 /*
2  * LALInferencePrior.c: Nested Sampling using LALInference
3  *
4  * Copyright (C) 2009 Ilya Mandel, Vivien Raymond, Christian Roever, Marc van der Sluys and 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 <lal/LALInferencePrior.h>
24 #include <lal/LALInferenceLikelihood.h>
25 #include <math.h>
26 #include <gsl/gsl_integration.h>
27 #include <gsl/gsl_cdf.h>
28 #include <gsl/gsl_errno.h>
29 #include <gsl/gsl_roots.h>
30 #include <lal/LALSimBurst.h>
31 #include <lal/XLALGSL.h>
32 
33 #include "logaddexp.h"
34 
35 #ifdef __GNUC__
36 #define UNUSED __attribute__ ((unused))
37 #else
38 #define UNUSED
39 #endif
40 
41 /* Private helper function prototypes */
42 static double qInnerIntegrand(double M2, void *viData);
43 static double etaInnerIntegrand(double M2, void *viData);
44 static double outerIntegrand(double M1, void *voData);
45 
46 static REAL8 REAL8max(REAL8 a, REAL8 b);
47 static REAL8 REAL8max(REAL8 a, REAL8 b)
48 {
49  return (a>b?a:b);
50 }
51 
52 void LALInferenceInitCBCPrior(LALInferenceRunState *runState)
53 {
54  char help[]="\
55  ----------------------------------------------\n\
56  --- Prior Arguments --------------------------\n\
57  ----------------------------------------------\n\
58  (--distance-prior-uniform) Impose uniform prior on distance and not volume (False)\n\
59  (--distance-prior-comoving-volume) Impose uniform prior on source-frame comoving volume (False)\n\
60  (--malmquistprior) Impose selection effects on the prior (False)\n\
61  (--malmquist-loudest-snr) Threshold SNR in the loudest detector (0.0)\n\
62  (--malmquist-second-loudest-snr) Threshold SNR in the second loudest detector (5.0)\n\
63  (--malmquist-network-snr) Threshold network SNR (0.0)\n\
64  (--analyticnullprior) Use analytic null prior\n\
65  (--nullprior) Use null prior in the sampled parameters\n\
66  (--alignedspin-zprior) Use prior on z component of spin that corresponds to fully precessing model\n\
67  (--spin-volumetricprior) Use prior on spin components that is uniform inside the sphere\n\
68  \n";
69  ProcessParamsTable *ppt = NULL;
70 
71  /* Print command line arguments if help requested */
72  if(runState == NULL || LALInferenceGetProcParamVal(runState->commandLine, "--help")) {
73  fprintf(stdout, "%s", help);
74  return;
75  }
76  ProcessParamsTable *commandLine=runState->commandLine;
77 
78  /* Choose the proper prior */
79  if (LALInferenceGetProcParamVal(commandLine, "--correlatedGaussianLikelihood") ||
80  LALInferenceGetProcParamVal(commandLine, "--bimodalGaussianLikelihood") ||
81  LALInferenceGetProcParamVal(commandLine, "--rosenbrockLikelihood") ||
82  LALInferenceGetProcParamVal(commandLine, "--analyticnullprior")) {
83  runState->prior = &LALInferenceAnalyticNullPrior;
84  runState->CubeToPrior = &LALInferenceAnalyticCubeToPrior;
85  } else if (LALInferenceGetProcParamVal(commandLine, "--nullprior")) {
86  runState->prior = &LALInferenceNullPrior;
87  /* CubeToPrior missing for null prior */
88  } else {
89  runState->prior = &LALInferenceInspiralPrior;
90  runState->CubeToPrior = &LALInferenceInspiralCubeToPrior;
91 
92  }
93 
94  /* Optional uniform prior on distance */
95  INT4 uniform_distance = 0;
96  INT4 src_comove_volume_distance = 0;
97  if (LALInferenceGetProcParamVal(commandLine, "--distance-prior-uniform")) {
98  uniform_distance = 1;
99  }
100  if (LALInferenceGetProcParamVal(commandLine, "--distance-prior-comoving-volume")) {
101  src_comove_volume_distance = 1;
102  }
103  LALInferenceAddVariable(runState->priorArgs,
104  "uniform_distance", &uniform_distance,
105  LALINFERENCE_INT4_t,
106  LALINFERENCE_PARAM_OUTPUT);
107 
108  LALInferenceAddVariable(runState->priorArgs,
109  "src_comove_volume_distance", &src_comove_volume_distance,
110  LALINFERENCE_INT4_t,
111  LALINFERENCE_PARAM_OUTPUT);
112 
113  /* Set up malmquist prior */
114  INT4 malmquist = 0;
115  if (LALInferenceGetProcParamVal(commandLine, "--malmquistprior")) {
116  malmquist = 1;
117  REAL8 malmquist_loudest = 0.0;
118  REAL8 malmquist_second_loudest = 5.0;
119  REAL8 malmquist_network = 0.0;
120 
121  ppt = LALInferenceGetProcParamVal(commandLine,
122  "--malmquist-loudest-snr");
123  if (ppt) malmquist_loudest = atof(ppt->value);
124 
125  ppt = LALInferenceGetProcParamVal(commandLine,
126  "--malmquist-second-loudest-snr");
127  if (ppt) malmquist_second_loudest = atof(ppt->value);
128 
129  ppt = LALInferenceGetProcParamVal(commandLine,
130  "--malmquist-network-snr");
131  if (ppt) malmquist_network = atof(ppt->value);
132 
133  LALInferenceAddVariable(runState->priorArgs,
134  "malmquist", &malmquist,
135  LALINFERENCE_INT4_t,
136  LALINFERENCE_PARAM_OUTPUT);
137 
138  LALInferenceAddVariable(runState->priorArgs,
139  "malmquist_loudest_snr",
140  &malmquist_loudest,
141  LALINFERENCE_REAL8_t,
142  LALINFERENCE_PARAM_OUTPUT);
143 
144  LALInferenceAddVariable(runState->priorArgs,
145  "malmquist_second_loudest_snr",
146  &malmquist_second_loudest,
147  LALINFERENCE_REAL8_t,
148  LALINFERENCE_PARAM_OUTPUT);
149 
150  LALInferenceAddVariable(runState->priorArgs,
151  "malmquist_network_snr",
152  &malmquist_network,
153  LALINFERENCE_REAL8_t,
154  LALINFERENCE_PARAM_OUTPUT);
155  }
156  INT4 one=1;
157  if(LALInferenceGetProcParamVal(commandLine,"--alignedspin-zprior"))
158  {
159  LALInferenceAddVariable(runState->priorArgs,"projected_aligned_spin",&one,LALINFERENCE_INT4_t,LALINFERENCE_PARAM_FIXED);
160  }
161  if(LALInferenceGetProcParamVal(commandLine,"--spin-volumetricprior"))
162  {
163  LALInferenceAddVariable(runState->priorArgs,"volumetric_spin",&one,LALINFERENCE_INT4_t,LALINFERENCE_PARAM_FIXED);
164  }
165  if(LALInferenceGetProcParamVal(commandLine,"--alignedspin-zprior")&&LALInferenceGetProcParamVal(commandLine,"--spin-volumetricprior"))
166  {
167  fprintf(stderr,"Error: You cannot use both --alignedspin-zprior and --spin-volumetricprior\n");
168  exit(1);
169  }
170  if(LALInferenceGetProcParamVal(commandLine,"--distance-prior-comoving-volume")&&LALInferenceGetProcParamVal(commandLine,"--distance-prior-uniform"))
171  {
172  fprintf(stderr,"Error: You cannot use both --distance-prior-comoving-volume and --distance-prior-uniform\n");
173  exit(1);
174  }
175 }
176 
177 void LALInferenceInitLIBPrior(LALInferenceRunState *runState)
178 {
179  /*Call CBC prior first, in case CBC approx is used, then check for burst approx and eventually overwrite runState->prior */
180  LALInferenceInitCBCPrior(runState);
181  /*LIB specific call in case of burst approximant */
182  ProcessParamsTable *commandLine=runState->commandLine;
183  ProcessParamsTable *ppt = NULL;
184  if ((ppt=LALInferenceGetProcParamVal(commandLine,"--approx"))){
185  if ((XLALCheckBurstApproximantFromString(ppt->value)))
186  {
187  /* Choose the proper prior */
188  if (LALInferenceGetProcParamVal(commandLine, "--correlatedGaussianLikelihood") ||
189  LALInferenceGetProcParamVal(commandLine, "--bimodalGaussianLikelihood") ||
190  LALInferenceGetProcParamVal(commandLine, "--analyticnullprior")) {
191  runState->prior = &LALInferenceAnalyticNullPrior;
192  } else if (LALInferenceGetProcParamVal(commandLine, "--nullprior")) {
193  runState->prior = &LALInferenceNullPrior;
194  } else {
195  runState->prior = &LALInferenceSineGaussianPrior;
196  }
197  }
198  }
199 }
200 
201 static REAL8 LALInferenceConstantCalibrationPrior(LALInferenceRunState *runState, LALInferenceVariables *params) {
202 
203  LALInferenceIFOData *ifo = NULL;
204  REAL8 ampWidth = -1.0;
205  REAL8 phaseWidth = -1.0;
206  REAL8 logPrior = 0.0;
207 
208  if (runState->commandLine == NULL || (!LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalAmp") &&
209  !LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalPha")))
210  {
211  return logPrior;
212  }
213 
214  ifo = runState->data;
215  do {
216 
217  if((LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalAmp"))){
218  ampWidth = *(REAL8 *)LALInferenceGetVariable(runState->priorArgs, "constcal_amp_uncertainty");
219  if (ampWidth>0){
220  char ampVarName[VARNAME_MAX];
221  REAL8 amp = 0.0;
222  if((VARNAME_MAX <= snprintf(ampVarName, VARNAME_MAX, "calamp_%s", ifo->name)))
223  {
224  fprintf(stderr,"variable name too long\n"); exit(1);
225  }
226  amp = *(REAL8*)LALInferenceGetVariable(params, ampVarName);
227  logPrior += -0.5*log(2.0*M_PI) - log(ampWidth) - 0.5*amp*amp/ampWidth/ampWidth;
228  }
229  }
230  if((LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalPha"))){
231  phaseWidth = *(REAL8 *)LALInferenceGetVariable(runState->priorArgs, "constcal_phase_uncertainty");
232  if (phaseWidth>0){
233  char phaseVarName[VARNAME_MAX];
234  REAL8 phase = 0.0;
235  if((VARNAME_MAX <= snprintf(phaseVarName, VARNAME_MAX, "calpha_%s", ifo->name)))
236  {
237  fprintf(stderr,"variable name too long\n"); exit(1);
238  }
239  phase = *(REAL8 *)LALInferenceGetVariable(params, phaseVarName);
240  logPrior += -0.5*log(2.0*M_PI) - log(phaseWidth) - 0.5*phase*phase/phaseWidth/phaseWidth;
241  }
242  }
243 
244  ifo = ifo->next;
245  } while (ifo);
246 
247  return logPrior;
248 }
249 
250 UINT4 LALInferenceCubeToConstantCalibrationPrior(LALInferenceRunState *runState, LALInferenceVariables *params, INT4 *idx, double *Cube, void UNUSED *context)
251 {
252  LALInferenceIFOData *ifo = NULL;
253  REAL8 ampWidth = -1.0;
254  REAL8 phaseWidth = -1.0;
255 
256  if (runState->commandLine == NULL || (!LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalAmp") &&
257  !LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalPha")))
258  {
259  return 1;
260  }
261  ifo = runState->data;
262  do {
263 
264  if((LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalAmp")))
265  {
266  ampWidth = *(REAL8 *)LALInferenceGetVariable(runState->priorArgs, "constcal_amp_uncertainty");
267  if (ampWidth>0)
268  {
269  char ampVarName[VARNAME_MAX];
270  REAL8 amp = 0.0;
271  if((VARNAME_MAX <= snprintf(ampVarName, VARNAME_MAX, "calamp_%s", ifo->name)))
272  {
273  fprintf(stderr,"variable name too long\n"); exit(1);
274  }
275  amp = LALInferenceCubeToGaussianPrior(Cube[(*idx)++], 0.0, ampWidth);
276  LALInferenceSetVariable(params, ampVarName, &amp);
277  }
278  }
279  if((LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalPha")))
280  {
281  phaseWidth = *(REAL8 *)LALInferenceGetVariable(runState->priorArgs, "constcal_phase_uncertainty");
282  if (phaseWidth>0)
283  {
284  char phaseVarName[VARNAME_MAX];
285  REAL8 phase = 0.0;
286  if((VARNAME_MAX <= snprintf(phaseVarName, VARNAME_MAX, "calpha_%s", ifo->name)))
287  {
288  fprintf(stderr,"variable name too long\n"); exit(1);
289  }
290  phase = LALInferenceCubeToGaussianPrior(Cube[(*idx)++], 0.0, phaseWidth);
291  LALInferenceSetVariable(params, phaseVarName, &phase);
292  }
293  }
294 
295  ifo = ifo->next;
296  } while (ifo);
297 
298  return 1;
299 
300 }
301 
302 
303 
304 /* Return the log Prior for the glitch amplitude */
305 REAL8 logGlitchAmplitudeDensity(REAL8 A, REAL8 Q, REAL8 f)
306 {
307  REAL8 SNR;
308  /*
309  REAL8 PIterm = 0.5*LAL_2_SQRTPI*LAL_SQRT1_2;
310  */
311  REAL8 PIterm = LAL_2_SQRTPI*LAL_SQRT1_2;
312  REAL8 SNRPEAK = 5.0;
313 
314  SNR = A*sqrt( (PIterm*Q/f) );
315  return log(SNR/(SNRPEAK*SNRPEAK))+(-SNR/SNRPEAK);
316 }
317 
318 /* Return the log Prior for the glitch model */
319 static REAL8 LALInferenceGlitchPrior(LALInferenceRunState *runState, LALInferenceVariables *params) {
320  LALInferenceVariableItem *item=params->head;
321  LALInferenceVariables *priorParams=runState->priorArgs;
322  REAL8 logPrior=0.0;
323 
324  /* check if glitch model is being used */
325  UINT4 glitchFlag = 0;
326  if(LALInferenceCheckVariable(params,"glitchFitFlag"))
327  glitchFlag = *((INT4 *)LALInferenceGetVariable(params, "glitchFitFlag"));
328 
329  if(glitchFlag)
330  {
331  UINT4 nifo,nglitch;
332  // Get parameters for current glitch model
333  UINT4Vector *gsize = *(UINT4Vector **) LALInferenceGetVariable(params, "glitch_size");
334 
335  gsl_matrix *glitch_f = *(gsl_matrix **)LALInferenceGetVariable(params, "morlet_f0");
336  gsl_matrix *glitch_Q = *(gsl_matrix **)LALInferenceGetVariable(params, "morlet_Q");
337  gsl_matrix *glitch_A = *(gsl_matrix **)LALInferenceGetVariable(params, "morlet_Amp");
338  gsl_matrix *gparams = NULL;
339 
340  REAL8 component_min=0.0;
341  REAL8 component_max=0.0;
342  REAL8 val=0.0;
343 
344  char priormin[512];
345  char priormax[512];
346 
347  REAL8 Anorm = *(REAL8 *)LALInferenceGetVariable(priorParams,"glitch_norm");
348 
349  REAL8 A,f,Q;
350  for(nifo=0; nifo<(UINT4)gsize->length; nifo++)
351  {
352  for(nglitch=0; nglitch<gsize->data[nifo]; nglitch++)
353  {
354  A = gsl_matrix_get(glitch_A,nifo,nglitch);
355  Q = gsl_matrix_get(glitch_Q,nifo,nglitch);
356  f = gsl_matrix_get(glitch_f,nifo,nglitch);
357 
358  logPrior += logGlitchAmplitudeDensity(A*Anorm,Q,f);
359  //printf("logPrior=%g\n",logPrior);
360  }
361  }
362  for(;item;item=item->next){
363  if(!strcmp(item->name,"morlet_f0" ) || !strcmp(item->name,"morlet_Q" ) || !strcmp(item->name,"morlet_t0" ) || !strcmp(item->name,"morlet_phi") )
364  {
365  gparams = *((gsl_matrix **)(item->value));
366 
367  sprintf(priormin,"%s_prior_min",item->name);
368  sprintf(priormax,"%s_prior_max",item->name);
369  component_min=*(REAL8 *)LALInferenceGetVariable(priorParams,priormin);
370  component_max=*(REAL8 *)LALInferenceGetVariable(priorParams,priormax);
371 
372  for(UINT4 i=0; i<gsize->length; i++)
373  {
374  for(UINT4 j=0; j<gsize->data[i]; j++)
375  {
376  val = gsl_matrix_get(gparams,i,j);
377 
378  //rejection sample on prior
379  if(val<component_min || val>component_max) return -INFINITY;
380  else logPrior -= log(component_max-component_min);
381  }
382  }
383  }//end morlet parameters prior
384  }//end loop through params
385  }//end check for glitch parameters
386 
387  return logPrior;
388 }
389 
390 /* Return the log Prior for the psd model */
391 static REAL8 LALInferencePSDPrior(LALInferenceRunState *runState, LALInferenceVariables *params)
392 {
393  LALInferenceVariables *priorParams=runState->priorArgs;
394  REAL8 logPrior=0.0;
395 
396 
397  /* check if PSD model is being used */
398  UINT4 psdFlag = 0;
399  if(LALInferenceCheckVariable(params, "psdScaleFlag"))
400  psdFlag = *((INT4 *)LALInferenceGetVariable(params, "psdScaleFlag"));
401 
402  /* PSD scale parameters */
403  if(psdFlag)
404  {
405  UINT4 i;
406  UINT4 j;
407 
408  REAL8 val = 0.0;
409  REAL8 var = 1.0;
410  REAL8 mean = 1.0;
411 
412  REAL8Vector *sigma = *((REAL8Vector **)LALInferenceGetVariable(priorParams, "psdsigma"));
413  gsl_matrix *nparams = *((gsl_matrix **)LALInferenceGetVariable(params, "psdscale"));
414 
415  REAL8 component_min=*(REAL8 *)LALInferenceGetVariable(priorParams,"psdrange_min");
416  REAL8 component_max=*(REAL8 *)LALInferenceGetVariable(priorParams,"psdrange_max");
417 
418  UINT4 psdGaussianPrior = *(UINT4 *)LALInferenceGetVariable(priorParams,"psdGaussianPrior");
419 
420  //Loop over IFOs
421  for(i=0; i<(UINT4)nparams->size1; i++)
422  {
423  //Loop over PSD windows
424  for(j=0; j<(UINT4)nparams->size2; j++)
425  {
426  var = sigma->data[j]*sigma->data[j];
427  val = gsl_matrix_get(nparams,i,j);
428 
429  //reject prior
430  if(val < component_min || val > component_max) return -INFINITY;
431  else if(psdGaussianPrior) logPrior += -0.5*( (mean-val)*(mean-val)/var + log(2.0*LAL_PI*var) );
432  }//end loop over windows
433  }//end loop over IFOs
434 
435  }//end psdFlag conditional
436 
437  return logPrior;
438 }
439 
440 
441 /* Return the log Prior of the variables specified, for the non-spinning/spinning inspiral signal case */
442 REAL8 LALInferenceInspiralPrior(LALInferenceRunState *runState, LALInferenceVariables *params, LALInferenceModel *model)
443 {
444  if (runState == NULL || runState->priorArgs == NULL || params == NULL)
445  XLAL_ERROR_REAL8(XLAL_EFAULT, "Null arguments received.");
446 
447  REAL8 logPrior=0.0;
448 
449  LALInferenceVariableItem *item=NULL;
450  LALInferenceVariables *priorParams=runState->priorArgs;
451  REAL8 min=-INFINITY, max=INFINITY;
452  REAL8 mc=0.0;
453  REAL8 m1=0.0,m2=0.0,q=0.0,eta=0.0;
454  REAL8 c0=1.012306, c1=1.136740, c2=0.262462, c3=0.016732, c4=0.000387; /* fitting coefficients for Will's cosmological distance prior, see https://git.ligo.org/RatesAndPopulations/lalinfsamplereweighting/blob/master/ApproxPrior.ipynb */
455 
456  /* check if signal model is being used */
457  UINT4 signalFlag=1;
458  if(LALInferenceCheckVariable(params, "signalModelFlag"))
459  signalFlag = *((INT4 *)LALInferenceGetVariable(params, "signalModelFlag"));
460 
461  if(signalFlag){
462 
463  /* Check boundaries for signal model parameters */
464  for(item=params->head;item;item=item->next)
465  {
466  if(item->vary==LALINFERENCE_PARAM_FIXED || item->vary==LALINFERENCE_PARAM_OUTPUT)
467  continue;
468  else if (LALInferenceCheckMinMaxPrior(priorParams, item->name))
469  {
470  if(item->type==LALINFERENCE_REAL8_t){
471  LALInferenceGetMinMaxPrior(priorParams, item->name, &min, &max);
472  if(*(REAL8 *) item->value < min || *(REAL8 *)item->value > max) return -INFINITY;
473  }
474  }
475  else if (LALInferenceCheckGaussianPrior(priorParams, item->name))
476  {
477  if(item->type==LALINFERENCE_REAL8_t){
478  REAL8 mean,stdev,val;
479  val = *(REAL8 *)item->value;
480  LALInferenceGetGaussianPrior(priorParams, item->name, &mean, &stdev);
481  logPrior+= -0.5*(mean-val)*(mean-val)/stdev/stdev - 0.5*log(LAL_TWOPI) - log(stdev);
482  }
483  }
484  }
485  if(LALInferenceCheckVariable(params, "flow") &&
486  LALInferenceCheckVariableNonFixed(params, "flow")) {
487  logPrior+=log(*(REAL8 *)LALInferenceGetVariable(params,"flow"));
488  }
489 
490 
491  if(LALInferenceCheckVariable(params,"logdistance"))
492  {
493  REAL8 log_dist = *(REAL8 *)LALInferenceGetVariable(params,"logdistance");
494  if ((LALInferenceCheckVariable(priorParams,"uniform_distance") && LALInferenceGetINT4Variable(priorParams,"uniform_distance"))) {
495  logPrior+=log_dist;
496  }
497  else if ((LALInferenceCheckVariable(priorParams,"src_comove_volume_distance") && LALInferenceGetINT4Variable(priorParams,"src_comove_volume_distance"))) {
498  REAL8 dist_Gpc = exp(log_dist)/1000.0;
499  REAL8 dist_Gpc2= dist_Gpc*dist_Gpc;
500  REAL8 dist_Gpc3= dist_Gpc2*dist_Gpc;
501  REAL8 dist_Gpc4= dist_Gpc3*dist_Gpc;
502  REAL8 denominator = c0+c1*dist_Gpc+c2*dist_Gpc2+c3*dist_Gpc3+c4*dist_Gpc4;
503  logPrior+=3.0* log_dist-log(denominator);
504  }
505  else {
506  logPrior+=3.0* log_dist;
507  }
508  }
509  else if(LALInferenceCheckVariable(params,"distance"))
510  {
511  if (!(LALInferenceCheckVariable(priorParams,"uniform_distance")&&LALInferenceGetINT4Variable(priorParams,"uniform_distance"))) {
512  REAL8 dist = *(REAL8 *)LALInferenceGetVariable(params,"distance");
513  if ((LALInferenceCheckVariable(priorParams,"src_comove_volume_distance") && LALInferenceGetINT4Variable(priorParams,"src_comove_volume_distance"))) {
514  REAL8 dist_Gpc = dist/1000.0;
515  REAL8 dist_Gpc2= dist_Gpc*dist_Gpc;
516  REAL8 dist_Gpc3= dist_Gpc2*dist_Gpc;
517  REAL8 dist_Gpc4= dist_Gpc3*dist_Gpc;
518  REAL8 denominator = c0+c1*dist_Gpc+c2*dist_Gpc2+c3*dist_Gpc3+c4*dist_Gpc4;
519  logPrior+=2.0*log(dist)-log(denominator);
520  }
521  else {
522  logPrior+=2.0*log(dist);
523  }
524  }
525  }
526  if(LALInferenceCheckVariable(params,"declination"))
527  {
528  /* Check that this is not an output variable */
529  if(LALInferenceGetVariableVaryType(params,"declination")==LALINFERENCE_PARAM_LINEAR)
530  logPrior+=log(fabs(cos(*(REAL8 *)LALInferenceGetVariable(params,"declination"))));
531  }
532 
533 
534  if(LALInferenceCheckVariable(params,"logmc")) {
535  mc=exp(*(REAL8 *)LALInferenceGetVariable(params,"logmc"));
536  } else if(LALInferenceCheckVariable(params,"chirpmass")) {
537  mc=(*(REAL8 *)LALInferenceGetVariable(params,"chirpmass"));
538  }
539 
540  if(LALInferenceCheckVariable(params,"q")) {
541  q=*(REAL8 *)LALInferenceGetVariable(params,"q");
542  LALInferenceMcQ2Masses(mc,q,&m1,&m2);
543  } else if(LALInferenceCheckVariable(params,"eta")) {
544  eta=*(REAL8 *)LALInferenceGetVariable(params,"eta");
545  LALInferenceMcEta2Masses(mc,eta,&m1,&m2);
546  }
547 
548  if(LALInferenceCheckVariable(params,"logmc")) {
549  if(LALInferenceCheckVariable(params,"q"))
550  logPrior+=log(m1*m1);
551  else
552  logPrior+=log(((m1+m2)*(m1+m2)*(m1+m2))/(m1-m2));
553  } else if(LALInferenceCheckVariable(params,"chirpmass")) {
554  if(LALInferenceCheckVariable(params,"q"))
555  logPrior+=log(m1*m1/mc);
556  else
557  logPrior+=log(((m1+m2)*(m1+m2))/((m1-m2)*pow(eta,3.0/5.0)));
558  }
559 
560  /* Check for individual mass priors */
561  if(LALInferenceCheckVariable(priorParams,"mass1_min"))
562  if(LALInferenceGetREAL8Variable(priorParams,"mass1_min") > m1)
563  return -INFINITY;
564  if(LALInferenceCheckVariable(priorParams,"mass1_max"))
565  if(LALInferenceGetREAL8Variable(priorParams,"mass1_max") < m1)
566  return -INFINITY;
567  if(LALInferenceCheckVariable(priorParams,"mass2_min"))
568  if(LALInferenceGetREAL8Variable(priorParams,"mass2_min") > m2)
569  return -INFINITY;
570  if(LALInferenceCheckVariable(priorParams,"mass2_max"))
571  if(LALInferenceGetREAL8Variable(priorParams,"mass2_max") < m2)
572  return -INFINITY;
573 
574 
575  if(LALInferenceCheckVariable(priorParams,"MTotMax"))
576  if(*(REAL8 *)LALInferenceGetVariable(priorParams,"MTotMax") < m1+m2)
577  return -INFINITY;
578 
579  if(LALInferenceCheckVariable(priorParams,"MTotMin"))
580  if(*(REAL8 *)LALInferenceGetVariable(priorParams,"MTotMin") > m1+m2)
581  return -INFINITY;
582 
583  if(model != NULL &&
584  LALInferenceCheckVariable(priorParams,"malmquist") &&
585  *(UINT4 *)LALInferenceGetVariable(priorParams,"malmquist") &&
586  !within_malmquist(runState, params, model))
587  return -INFINITY;
588 
589  UINT4 volumetric_spins = LALInferenceCheckVariable(runState->priorArgs,"volumetric_spin") && LALInferenceGetVariable(runState->priorArgs,"volumetric_spin");
590  /* Apply spin priors for precessing case */
591  if(LALInferenceCheckVariable(params,"tilt_spin1"))
592  {
593  LALInferenceParamVaryType vtype=LALInferenceGetVariableVaryType(params,"tilt_spin1");
594  if(vtype!=LALINFERENCE_PARAM_FIXED && vtype!=LALINFERENCE_PARAM_OUTPUT)
595  {
596  if(volumetric_spins)
597  {
598  /* homogenous inside spin bound */
599  /* V = (4/3)*pi*(a_max^3 - a_min^3) */
600  REAL8 a = LALInferenceGetREAL8Variable(params,"a_spin1");
601  REAL8 a_max,a_min;
602  LALInferenceGetMinMaxPrior(runState->priorArgs,"a_spin1",&a_min,&a_max);
603  REAL8 V = (4./3.)*LAL_PI * (a_max*a_max*a_max - a_min*a_min*a_min);
604  logPrior+=log(fabs(a*a))-log(fabs(V));
605  }
606  /* Usual case has uniform in a, but both cases have sin(tilt) from volume element */
607  logPrior+=log(fabs(sin(*(REAL8 *)LALInferenceGetVariable(params,"tilt_spin1"))));
608  }
609  }
610  else
611  {
612  if(volumetric_spins)
613  {
614  /* Volumetric prior marginalised onto z component */
615  REAL8 a = LALInferenceGetREAL8Variable(params,"a_spin1");
616  REAL8 a_max,a_min;
617  LALInferenceGetMinMaxPrior(runState->priorArgs,"a_spin1",&a_min,&a_max);
618  REAL8 V = (4./3.)*LAL_PI * (a_max*a_max*a_max);
619  logPrior+=log(fabs((3./4.)*(a_max*a_max - a*a)))-log(fabs(V));
620  }
621  }
622  if(LALInferenceCheckVariable(params,"tilt_spin2"))
623  {
624  LALInferenceParamVaryType vtype=LALInferenceGetVariableVaryType(params,"tilt_spin2");
625  if(vtype!=LALINFERENCE_PARAM_FIXED && vtype!=LALINFERENCE_PARAM_OUTPUT)
626  {
627  if(volumetric_spins)
628  {
629  REAL8 a = LALInferenceGetREAL8Variable(params,"a_spin2");
630  REAL8 a_max,a_min;
631  LALInferenceGetMinMaxPrior(runState->priorArgs,"a_spin2",&a_min,&a_max);
632  REAL8 V = (4./3.)*LAL_PI * (a_max*a_max*a_max - a_min*a_min*a_min);
633  logPrior+=log(fabs(a*a))-log(fabs(V));
634  }
635  logPrior+=log(fabs(sin(*(REAL8 *)LALInferenceGetVariable(params,"tilt_spin2"))));
636  }
637  }
638  else
639  {
640  if(volumetric_spins)
641  {
642  /* Volumetric prior marginalised onto z component */
643  REAL8 a = LALInferenceGetREAL8Variable(params,"a_spin2");
644  REAL8 a_max,a_min;
645  LALInferenceGetMinMaxPrior(runState->priorArgs,"a_spin2",&a_min,&a_max);
646  REAL8 V = (4./3.)*LAL_PI * (a_max*a_max*a_max);
647  logPrior+=log(fabs((3./4.)*(a_max*a_max - a*a)))-log(fabs(V));
648  }
649  }
650 
651  /* Optional prior on aligned spin component that corresponds to the effective prior on
652  that component when using a precessing spin model. p(z) = (1/2)(1/R)log(|z|/R)
653  Where R is the maximum magnitude of the spin vector max(|a_spin1_max|,|a_spin1_min|).
654  */
655  if (LALInferenceCheckVariable(priorParams,"projected_aligned_spin") && LALInferenceGetVariable(priorParams,"projected_aligned_spin"))
656  {
657  REAL8 z=0.0;
658  /* Double-check for tilts to prevent accidental double-prior */
659  if(LALInferenceCheckVariable(params,"a_spin1") && !LALInferenceCheckVariable(params,"tilt_spin1"))
660  {
661  REAL8 R = REAL8max(fabs(LALInferenceGetREAL8Variable(priorParams,"a_spin1_max")),fabs(LALInferenceGetREAL8Variable(priorParams,"a_spin1_min")));
662  z=LALInferenceGetREAL8Variable(params,"a_spin1");
663  logPrior += -log(2.0) - log(R) + log(-log(fabs(z) / R));
664  }
665  if(LALInferenceCheckVariable(params,"a_spin2")&& !LALInferenceCheckVariable(params,"tilt_spin2"))
666  {
667  REAL8 R = REAL8max(fabs(LALInferenceGetREAL8Variable(priorParams,"a_spin2_max")),fabs(LALInferenceGetREAL8Variable(priorParams,"a_spin2_min")));
668  z=LALInferenceGetREAL8Variable(params,"a_spin2");
669  logPrior += -log(2.0) - log(R) + log(-log(fabs(z) / R));
670  }
671 
672  }
673 
674  // Apply additional prior if using eos parameters
675  if((LALInferenceCheckVariable(params,"logp1")&&LALInferenceCheckVariable(params,"gamma1")&&LALInferenceCheckVariable(params,"gamma2")&&LALInferenceCheckVariable(params,"gamma3")))
676  {
677  /*If EOS params and masses are aphysical, return -INFINITY to ensure point is rejected*/
678  if(LALInferenceEOSPhysicalCheck(params,runState->commandLine)==XLAL_FAILURE){
679  return -INFINITY;
680  }
681  }
682  else if((LALInferenceCheckVariable(params,"SDgamma0")&&LALInferenceCheckVariable(params,"SDgamma1")&&LALInferenceCheckVariable(params,"SDgamma2")&&LALInferenceCheckVariable(params,"SDgamma3")))
683  {
684  /*If EOS params and masses are aphysical, return -INFINITY to ensure point is rejected*/
685  if(LALInferenceEOSPhysicalCheck(params,runState->commandLine)==XLAL_FAILURE){
686  return -INFINITY;
687  }
688  }
689 
690  }/* end prior for signal model parameters */
691 
692 
693  /* Calibration priors. */
694  /* Disabled as this is now handled automatically */
695  //logPrior += LALInferenceSplineCalibrationPrior(runState, params);
696  logPrior += LALInferenceConstantCalibrationPrior(runState, params);
697  /* Evaluate PSD prior (returns 0 if no PSD model) */
698  logPrior += LALInferencePSDPrior(runState, params);
699 
700  /* Evaluate glitch prior (returns 0 if no glitch model) */
701  logPrior += LALInferenceGlitchPrior(runState, params);
702 
703  return(logPrior);
704 }
705 
706 /* Convert the hypercube parameter to physical parameters, for the non-spinning inspiral signal case */
707 UINT4 LALInferenceInspiralCubeToPrior(LALInferenceRunState *runState, LALInferenceVariables *params, LALInferenceModel *model, double *Cube, void *context)
708 {
709  REAL8 min=-INFINITY, max=INFINITY;
710  LALInferenceVariableItem *item;
711  LALInferenceVariables *priorParams=runState->priorArgs;
712 
713  char **info = (char **)context;
714  char *timeID = &info[2][0];
715  int i = 0;
716 
717  INT4 SKY_FRAME=0;
718  if(LALInferenceCheckVariable(params,"SKY_FRAME"))
719  SKY_FRAME=*(INT4 *)LALInferenceGetVariable(params,"SKY_FRAME");
720  double azimuth, cosalpha, lat, longitude, t0, tc;
721 
722  if(SKY_FRAME==1)
723  {
724  // detector frame azimuth
725  if (LALInferenceCheckVariable(params, "azimuth"))
726  {
727  item = LALInferenceGetItem(params, "azimuth");
728  if (item->vary != LALINFERENCE_PARAM_FIXED)
729  {
730  LALInferenceGetMinMaxPrior(runState->priorArgs, "azimuth", (void *)&min, (void *)&max);
731  azimuth = LALInferenceCubeToFlatPrior(Cube[i], min, max);
732  LALInferenceSetVariable(params, "azimuth", &azimuth);
733  i++;
734  }
735  }
736 
737  // detector frame cosalpha
738  if (LALInferenceCheckVariable(params, "cosalpha"))
739  {
740  item = LALInferenceGetItem(params, "cosalpha");
741  if (item->vary != LALINFERENCE_PARAM_FIXED)
742  {
743  LALInferenceGetMinMaxPrior(runState->priorArgs, "cosalpha", (void *)&min, (void *)&max);
744  cosalpha = LALInferenceCubeToFlatPrior(Cube[i], min, max);
745  LALInferenceSetVariable(params, "cosalpha", &cosalpha);
746  i++;
747  }
748  }
749  }
750  else
751  {
752  // latitude
753  if (LALInferenceCheckVariable(params, "declination"))
754  {
755  item = LALInferenceGetItem(params, "declination");
756  if(item->vary != LALINFERENCE_PARAM_FIXED)
757  {
758  lat = asin(2.0 * Cube[i] - 1.0);
759  LALInferenceSetVariable(params, "declination", &lat);
760  i++;
761  }
762  }
763 
764  // longitude
765  if (LALInferenceCheckVariable(params, "rightascension"))
766  {
767  item = LALInferenceGetItem(params, "rightascension");
768  if(item->vary != LALINFERENCE_PARAM_FIXED)
769  {
770  LALInferenceGetMinMaxPrior(runState->priorArgs, "rightascension", (void *)&min, (void *)&max);
771  longitude = LALInferenceCubeToFlatPrior(Cube[i], min, max);
772  LALInferenceSetVariable(params, "rightascension", &longitude);
773  i++;
774  }
775  }
776  }
777 
778  // phi
779  if (LALInferenceCheckVariable(params, "phase"))
780  {
781  item = LALInferenceGetItem(params, "phase");
782  if(item->vary != LALINFERENCE_PARAM_FIXED)
783  {
784  LALInferenceGetMinMaxPrior(runState->priorArgs, "phase", (void *)&min, (void *)&max);
785  double phi = LALInferenceCubeToFlatPrior(Cube[i], min, max);
786  LALInferenceSetVariable(params, "phase", &phi);
787  i++;
788  }
789  }
790 
791  // psi
792  if (LALInferenceCheckVariable(params, "polarisation"))
793  {
794  item = LALInferenceGetItem(params, "polarisation");
795  if(item->vary != LALINFERENCE_PARAM_FIXED)
796  {
797  LALInferenceGetMinMaxPrior(runState->priorArgs, "polarisation", (void *)&min, (void *)&max);
798  double psi = LALInferenceCubeToFlatPrior(Cube[i], min, max);
799  LALInferenceSetVariable(params, "polarisation", &psi);
800  i++;
801  }
802  }
803 
804  if(SKY_FRAME==1)
805  {
806  // detector frame t0
807  if (LALInferenceCheckVariable(params, "t0"))
808  {
809  item = LALInferenceGetItem(params, "t0");
810  if(item->vary != LALINFERENCE_PARAM_FIXED)
811  {
812  LALInferenceGetMinMaxPrior(runState->priorArgs, "t0", (void *)&min, (void *)&max);
813  t0 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
814  LALInferenceSetVariable(params, "t0", &t0);
815  sprintf(timeID,"%d",i);
816  i++;
817  }
818  }
819  }
820  else
821  {
822  // time
823  if (LALInferenceCheckVariable(params, "time"))
824  {
825  item = LALInferenceGetItem(params, "time");
826  if(item->vary != LALINFERENCE_PARAM_FIXED)
827  {
828  LALInferenceGetMinMaxPrior(runState->priorArgs, "time", (void *)&min, (void *)&max);
829  tc = LALInferenceCubeToFlatPrior(Cube[i], min, max);
830  LALInferenceSetVariable(params, "time", &tc);
831  sprintf(timeID,"%d",i);
832  i++;
833  }
834  }
835  }
836 
837 
838  // mass variables
839  double mc = 0.0, eta = 0.0, q = 0.0, m1 = 0.0, m2 = 0.0, m = 0.0;
840 
841  // check if mchirp is fixed
842  if( LALInferenceCheckVariable(params,"logmc") )
843  {
844  item = LALInferenceGetItem(params, "logmc");
845  if(item->vary == LALINFERENCE_PARAM_FIXED) mc = exp(*(REAL8 *)LALInferenceGetVariable(params, "logmc"));
846  }
847  else if( LALInferenceCheckVariable(params,"chirpmass") )
848  {
849  item = LALInferenceGetItem(params, "chirpmass");
850  if(item->vary == LALINFERENCE_PARAM_FIXED) mc = *(REAL8 *)LALInferenceGetVariable(params, "chirpmass");
851  }
852 
853  // check if eta is fixed
854  if( LALInferenceCheckVariable(params,"eta") )
855  {
856  item = LALInferenceGetItem(params, "eta");
857  if(item->vary == LALINFERENCE_PARAM_FIXED)
858  {
859  eta = *(REAL8 *)LALInferenceGetVariable(params, "eta");
860  if( mc != 0.0 ) LALInferenceMcEta2Masses(mc,eta,&m1,&m2);
861  }
862  }
863  else if( LALInferenceCheckVariable(params,"q") )
864  {
865  item = LALInferenceGetItem(params, "q");
866  if(item->vary == LALINFERENCE_PARAM_FIXED)
867  {
868  q = *(REAL8 *)LALInferenceGetVariable(params, "q");
869  if( mc != 0.0 ) LALInferenceMcQ2Masses(mc,q,&m1,&m2);
870  }
871  }
872 
873  //m1 & m2
874  if( m1 == 0.0 && m2 == 0.0 )
875  {
876  REAL8 m1_min = *(REAL8 *)LALInferenceGetVariable(priorParams,"mass1_min");
877  REAL8 m2_min = *(REAL8 *)LALInferenceGetVariable(priorParams,"mass2_min");
878  REAL8 m1_max = *(REAL8 *)LALInferenceGetVariable(priorParams,"mass1_max");
879  REAL8 m2_max = *(REAL8 *)LALInferenceGetVariable(priorParams,"mass2_max");
880 
881  if( m1_min == m1_max && m2_min==m2_max)
882  {
883  m1 = m1_min;
884  m2 = m2_min;
885  m = m1 + m2;
886  eta = m1 * m2 / (m*m);
887  mc = pow(eta,0.6) * m;
888  q = m2 / m1; // asymmetric mass ratio, m1 >= m2
889  }
890  else
891  {
892  m1 = LALInferenceCubeToFlatPrior(Cube[i], m1_min, m1_max);
893  m2 = LALInferenceCubeToFlatPrior(Cube[i+1], m2_min, m2_max);
894  if(m1<m2)
895  {
896  double temp = m2;
897  m2 = m1;
898  m1 = temp;
899  }
900 
901  m = m1 + m2;
902  eta = m1 * m2 / (m*m);
903  mc = pow(eta,0.6) * m;
904  q = m2 / m1; // asymmetric mass ratio, m1 >= m2
905  i++;
906  i++;
907  }
908 
909  // chirp mass and eta/q
910  if(LALInferenceCheckVariable(params,"eta")||LALInferenceCheckVariable(params,"q"))
911  {
912  if(LALInferenceCheckVariable(params,"logmc"))
913  {
914  double logmc = log(mc);
915  LALInferenceSetVariable(params, "logmc", &logmc);
916  }
917  else if(LALInferenceCheckVariable(params,"chirpmass"))
918  {
919  LALInferenceSetVariable(params, "chirpmass", &mc);
920  }
921 
922  if(LALInferenceCheckVariable(params,"q"))
923  LALInferenceSetVariable(params, "q", &q);
924  else if(LALInferenceCheckVariable(params,"eta"))
925  LALInferenceSetVariable(params, "eta", &eta);
926  }
927  }
928 
929  // distance
930  double dist;
931  if( LALInferenceCheckVariable(params,"logdistance") )
932  {
933  item = LALInferenceGetItem(params, "logdistance");
934  if(item->vary != LALINFERENCE_PARAM_FIXED)
935  {
936  LALInferenceGetMinMaxPrior(runState->priorArgs, "logdistance", (void *)&min, (void *)&max);
937  min = exp(min); max = exp(max);
938  dist = LALInferenceCubeToPowerPrior(2.0, Cube[i], min, max);
939  double logdist = log(dist);
940  LALInferenceSetVariable(params, "logdistance", &logdist);
941  i++;
942  }
943  }
944  else if( LALInferenceCheckVariable(params,"distance") )
945  {
946  item = LALInferenceGetItem(params, "distance");
947  if(item->vary != LALINFERENCE_PARAM_FIXED)
948  {
949  LALInferenceGetMinMaxPrior(runState->priorArgs, "distance", (void *)&min, (void *)&max);
950  dist = LALInferenceCubeToPowerPrior(2.0, Cube[i], min, max);
951  LALInferenceSetVariable(params, "distance", &dist);
952  i++;
953  }
954  }
955 
956  // a_spin1
957  if(LALInferenceCheckVariable(params,"a_spin1"))
958  {
959  item = LALInferenceGetItem(params, "a_spin1");
960  if(item->vary != LALINFERENCE_PARAM_FIXED)
961  {
962  LALInferenceGetMinMaxPrior(runState->priorArgs, "a_spin1", (void *)&min, (void *)&max);
963  double a_spin1 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
964  LALInferenceSetVariable(params, "a_spin1", &a_spin1);
965  i++;
966  }
967  }
968  // a_spin2
969  if(LALInferenceCheckVariable(params,"a_spin2"))
970  {
971  item = LALInferenceGetItem(params, "a_spin2");
972  if(item->vary != LALINFERENCE_PARAM_FIXED)
973  {
974  LALInferenceGetMinMaxPrior(runState->priorArgs, "a_spin2", (void *)&min, (void *)&max);
975  double a_spin2 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
976  LALInferenceSetVariable(params, "a_spin2", &a_spin2);
977  i++;
978  }
979  }
980 
981  // cos theta_JN for system-frame parameters
982  if(LALInferenceCheckVariable(params,"costheta_jn"))
983  {
984  item = LALInferenceGetItem(params, "costheta_jn");
985  if(item->vary != LALINFERENCE_PARAM_FIXED)
986  {
987  LALInferenceGetMinMaxPrior(runState->priorArgs, "costheta_jn", (void *)&min, (void *)&max);
988  double costheta_JN = LALInferenceCubeToFlatPrior(Cube[i], min, max);
989  LALInferenceSetVariable(params, "costheta_jn", &costheta_JN);
990  i++;
991  }
992  }
993 
994  // phi_JL for system-frame parameters
995  if(LALInferenceCheckVariable(params,"phi_jl"))
996  {
997  item = LALInferenceGetItem(params, "phi_jl");
998  if(item->vary != LALINFERENCE_PARAM_FIXED)
999  {
1000  LALInferenceGetMinMaxPrior(runState->priorArgs, "phi_jl", (void *)&min, (void *)&max);
1001  double phi_JL = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1002  LALInferenceSetVariable(params, "phi_jl", &phi_JL);
1003  i++;
1004  }
1005  }
1006 
1007  // tilt of spin 1 for system-frame parameters
1008  if(LALInferenceCheckVariable(params,"tilt_spin1"))
1009  {
1010  item = LALInferenceGetItem(params, "tilt_spin1");
1011  if(item->vary != LALINFERENCE_PARAM_FIXED)
1012  {
1013  LALInferenceGetMinMaxPrior(runState->priorArgs, "tilt_spin1", (void *)&min, (void *)&max);
1014  double tilt_spin1 = LALInferenceCubeToSinPrior(Cube[i], min, max);
1015  LALInferenceSetVariable(params, "tilt_spin1", &tilt_spin1);
1016  i++;
1017  }
1018  }
1019 
1020  // tilt of spin 2 for system-frame parameters
1021  if(LALInferenceCheckVariable(params,"tilt_spin2"))
1022  {
1023  item = LALInferenceGetItem(params, "tilt_spin2");
1024  if(item->vary != LALINFERENCE_PARAM_FIXED)
1025  {
1026  LALInferenceGetMinMaxPrior(runState->priorArgs, "tilt_spin2", (void *)&min, (void *)&max);
1027  double tilt_spin2 = LALInferenceCubeToSinPrior(Cube[i], min, max);
1028  LALInferenceSetVariable(params, "tilt_spin2", &tilt_spin2);
1029  i++;
1030  }
1031  }
1032 
1033  // phi12 for system-frame parameters
1034  if(LALInferenceCheckVariable(params,"phi12"))
1035  {
1036  item = LALInferenceGetItem(params, "phi12");
1037  if(item->vary != LALINFERENCE_PARAM_FIXED)
1038  {
1039  LALInferenceGetMinMaxPrior(runState->priorArgs, "phi12", (void *)&min, (void *)&max);
1040  double phi12 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1041  LALInferenceSetVariable(params, "phi12", &phi12);
1042  i++;
1043  }
1044  }
1045 
1046  UINT4 ScaleTest = LALInferenceCubeToPSDScaleParams(priorParams, params, &i, Cube, context);
1047  UINT4 ConstCalib = LALInferenceCubeToConstantCalibrationPrior(runState, params, &i, Cube, context);
1048  //UINT4 SplineCalib = LALInferenceCubeToSplineCalibrationPrior(runState, params, &i, Cube, context);
1049 
1050  /* Check boundaries */
1051  if (ScaleTest==0 || ConstCalib==0 /*|| SplineCalib==0*/) return 0;
1052  item=params->head;
1053  for(;item;item=item->next)
1054  {
1055  if(item->vary==LALINFERENCE_PARAM_FIXED || item->vary==LALINFERENCE_PARAM_OUTPUT)
1056  continue;
1057  else if (item->type != LALINFERENCE_gslMatrix_t && item->type != LALINFERENCE_REAL8Vector_t)
1058  {
1059  LALInferenceGetMinMaxPrior(priorParams, item->name, (void *)&min, (void *)&max);
1060  if(*(REAL8 *) item->value < min || *(REAL8 *)item->value > max) return 0 ;
1061  }
1062  }
1063 
1064  /* Check for real mass values */
1065  if(LALInferenceCheckVariable(params,"logmc"))
1066  if(isnan(*(REAL8 *)LALInferenceGetVariable(params,"logmc")))
1067  return 0;
1068  if(LALInferenceCheckVariable(params,"chirpmass"))
1069  if(isnan(*(REAL8 *)LALInferenceGetVariable(params,"chirpmass")))
1070  return 0;
1071  if(LALInferenceCheckVariable(params,"eta"))
1072  if(isnan(*(REAL8 *)LALInferenceGetVariable(params,"eta"))
1073  ||*(REAL8 *)LALInferenceGetVariable(params,"eta") < 0.0
1074  || *(REAL8 *)LALInferenceGetVariable(params,"eta") > 0.25)
1075  return 0;
1076  if(LALInferenceCheckVariable(params,"q"))
1077  if(isnan(*(REAL8 *)LALInferenceGetVariable(params,"q"))
1078  ||*(REAL8 *)LALInferenceGetVariable(params,"q") < 0.0
1079  || *(REAL8 *)LALInferenceGetVariable(params,"q") > 1.0)
1080  return 0;
1081 
1082  if(LALInferenceCheckVariable(priorParams,"MTotMin"))
1083  if(*(REAL8 *)LALInferenceGetVariable(priorParams,"MTotMin") > m1+m2)
1084  return 0;
1085 
1086  if(LALInferenceCheckVariable(priorParams,"MTotMax"))
1087  if(*(REAL8 *)LALInferenceGetVariable(priorParams,"MTotMax") < m1+m2)
1088  return 0;
1089 
1090  /* Check for individual mass priors */
1091  if(LALInferenceCheckVariable(priorParams,"mass1_min"))
1092  if(LALInferenceGetREAL8Variable(priorParams,"mass1_min") > m1)
1093  return 0;
1094  if(LALInferenceCheckVariable(priorParams,"mass1_max"))
1095  if(LALInferenceGetREAL8Variable(priorParams,"mass1_max") < m1)
1096  return 0;
1097  if(LALInferenceCheckVariable(priorParams,"mass2_min"))
1098  if(LALInferenceGetREAL8Variable(priorParams,"mass2_min") > m2)
1099  return 0;
1100  if(LALInferenceCheckVariable(priorParams,"mass2_max"))
1101  if(LALInferenceGetREAL8Variable(priorParams,"mass2_max") < m2)
1102  return 0;
1103 
1104  if(model != NULL &&
1105  LALInferenceCheckVariable(priorParams,"malmquist") &&
1106  *(UINT4 *)LALInferenceGetVariable(priorParams,"malmquist") &&
1107  !within_malmquist(runState, params, model))
1108  return 0;
1109 
1110  return 1;
1111 }
1112 
1113 void LALInferenceCyclicReflectiveBound(LALInferenceVariables *parameter,
1114  LALInferenceVariables *priorArgs){
1115  REAL8 val, offset, delta, min, max;
1116  UINT8 n;
1117 
1118  if (parameter == NULL || priorArgs == NULL)
1119  XLAL_ERROR_VOID(XLAL_EFAULT, "Null arguments received.");
1120 
1121  /* Apply cyclic and reflective boundaries to parameter to bring it back
1122  within the prior */
1123  LALInferenceVariableItem *paraHead=NULL;
1124 
1125  for (paraHead=parameter->head; paraHead; paraHead=paraHead->next) {
1126  if( paraHead->vary==LALINFERENCE_PARAM_FIXED ||
1127  paraHead->vary==LALINFERENCE_PARAM_OUTPUT ||
1128  !LALInferenceCheckMinMaxPrior(priorArgs, paraHead->name) ) continue;
1129 
1130  LALInferenceGetMinMaxPrior(priorArgs,paraHead->name, &min, &max);
1131 
1132  /* Check that the minimum and maximum make sense. */
1133  if (min >= max) {
1134  XLAL_ERROR_VOID(XLAL_EINVAL, "Minimum %f for variable '%s' is not less than maximum %f.", min, paraHead->name, max);
1135  }
1136  delta = max - min;
1137 
1138  val = *(REAL8 *)paraHead->value;
1139 
1140  if (val == INFINITY)
1141  return;
1142 
1143  // Nothing to do if between bounds
1144  if ((val >= min) && (val <= max))
1145  continue;
1146 
1147  if(paraHead->vary==LALINFERENCE_PARAM_CIRCULAR) {
1148  /* For cyclic boundaries, mod out by range. */
1149  if (val > max) {
1150  offset = val - min;
1151  *(REAL8 *)paraHead->value = min + fmod(offset, delta);
1152  } else {
1153  offset = max - val;
1154  *(REAL8 *)paraHead->value = max - fmod(offset, delta);
1155  }
1156  } else if (paraHead->vary==LALINFERENCE_PARAM_LINEAR && paraHead->type==LALINFERENCE_REAL8_t) {
1157  /* For linear boundaries, reflect about endpoints of range until
1158  within range. SKIP NOISE PARAMETERS (ONLY CHECK REAL8) */
1159  /*
1160  This reflection is accomplished by also modding out the 'excess' (that is,
1161  the difference between the value and the upper or lower bound, as appropriate)
1162  as done above for cyclic parameters. However, unlike in the cyclic case, to actually
1163  agree with what the reflection would give, we must also note whether the range divides
1164  into the excess an even or an odd number of times, and therefore whether we add the
1165  remainder to the lower bound or subtract it from the upper bound. Which we do also
1166  depends on whether the value was below the minimum or above the maximum.
1167  */
1168  if (val > max) {
1169  offset = val - max;
1170  // How many times do we fold?
1171  n = (UINT8) (offset/delta);
1172  // Now make offset the remainder
1173  offset = fmod(offset, delta);
1174  if (n % 2){
1175  // If we fold an odd number of times,
1176  // then add offset to min
1177  *(REAL8 *)paraHead->value = min + offset;
1178  } else {
1179  // Otherwise, subtract it from max
1180  *(REAL8 *)paraHead->value = max - offset;
1181  }
1182  } else {
1183  // We only get here if val < min
1184  offset = min - val;
1185  // How many times do we fold?
1186  n = (UINT8) (offset/delta);
1187  // Now make offset the remainder
1188  offset = fmod(offset, delta);
1189  if (n % 2){
1190  // If we fold an odd number of times,
1191  // then subtract offset from max
1192  *(REAL8 *)paraHead->value = max - offset;
1193  } else {
1194  // Otherwise, add it to min
1195  *(REAL8 *)paraHead->value = min + offset;
1196  }
1197  }
1198  }
1199  }
1200  return;
1201 }
1202 
1203 
1204 void LALInferenceRotateInitialPhase( LALInferenceVariables *parameter){
1205 
1206  if (parameter == NULL)
1207  XLAL_ERROR_VOID(XLAL_EFAULT, "Null arguments received.");
1208 
1209  LALInferenceVariableItem *paraHead = NULL;
1210  LALInferenceVariableItem *paraPhi0 = parameter->head;
1211  REAL8 rotphi0 = 0.;
1212  UINT4 idx1 = 0, idx2 = 0;
1213 
1214  for(paraHead=parameter->head;paraHead;paraHead=paraHead->next){
1215  if (paraHead->vary == LALINFERENCE_PARAM_CIRCULAR &&
1216  !strcmp(paraHead->name, "psi") ){
1217  /* if psi is outside the -pi/4 -> pi/4 boundary the set to rotate phi0
1218  by pi (psi will have been rescaled to be between 0 to 2pi as a
1219  circular parameter).*/
1220  if (*(REAL8 *)paraHead->value > LAL_TWOPI ||
1221  *(REAL8 *)paraHead->value < 0. ) rotphi0 = LAL_PI;
1222 
1223  /* Psi has been found; set flag. */
1224  idx1++;
1225  }
1226 
1227  /* If phi0 is found, set some flags to indicate so. */
1228  if (paraHead->vary == LALINFERENCE_PARAM_CIRCULAR &&
1229  !strcmp(paraHead->name, "phi0") ){
1230  idx1++;
1231  idx2++;
1232  }
1233 
1234  /* Advance paraPhi0 if phi0 hasn't yet been found. */
1235  if ( idx2 == 0 ) paraPhi0 = paraPhi0->next;
1236 
1237  /* If both variables have been found, stop here. */
1238  if ( idx1 == 2 ) break;
1239  }
1240 
1241  if( rotphi0 != 0. ) *(REAL8 *)paraPhi0->value += rotphi0;
1242 
1243  return;
1244 }
1245 
1246 /* Return the log Prior of the variables specified for the sky localisation project, ref: https://www.lsc-group.phys.uwm.edu/ligovirgo/cbcnote/SkyLocComparison#priors, for the non-spinning/spinning inspiral signal case */
1247 REAL8 LALInferenceInspiralSkyLocPrior(LALInferenceRunState *runState, LALInferenceVariables *params, UNUSED LALInferenceModel *model)
1248 {
1249  REAL8 logPrior=0.0;
1250  REAL8 val=0.0;
1251  static int SkyLocPriorWarning = 0;
1252  (void)runState;
1253  LALInferenceVariableItem *item=NULL;
1254  LALInferenceVariables *priorParams=runState->priorArgs;
1255  REAL8 min=-INFINITY, max=INFINITY;
1256  REAL8 logmc=0.0,mc=0.0;
1257  REAL8 m1=0.0,m2=0.0,q=0.0,eta=0.0;
1258 
1259  if (!SkyLocPriorWarning ) {
1260  SkyLocPriorWarning = 1;
1261  fprintf(stderr, "SkyLocalization priors are being used. (in %s, line %d)\n", __FILE__, __LINE__);
1262  }
1263  /* Check boundaries for signal model parameters */
1264  for(item=params->head;item;item=item->next)
1265  {
1266  if(item->vary==LALINFERENCE_PARAM_FIXED || item->vary==LALINFERENCE_PARAM_OUTPUT)
1267  continue;
1268  else if (LALInferenceCheckMinMaxPrior(priorParams, item->name))
1269  {
1270  if(item->type==LALINFERENCE_REAL8_t){
1271  LALInferenceGetMinMaxPrior(priorParams, item->name, &min, &max);
1272  if(*(REAL8 *) item->value < min || *(REAL8 *)item->value > max) return -INFINITY;
1273  }
1274  }
1275  else if (LALInferenceCheckGaussianPrior(priorParams, item->name))
1276  {
1277  if(item->type==LALINFERENCE_REAL8_t){
1278  REAL8 mean,stdev;
1279  val = *(REAL8 *)item->value;
1280  LALInferenceGetGaussianPrior(priorParams, item->name, &mean, &stdev);
1281  logPrior+= -0.5*(mean-val)*(mean-val)/stdev/stdev - 0.5*log(LAL_TWOPI) - log(stdev);
1282  }
1283  }
1284  }
1285 
1286  /*Use a uniform in log D distribution*/
1287 
1288  if(LALInferenceCheckVariable(params, "flow") &&
1289  LALInferenceCheckVariableNonFixed(params, "flow")) {
1290  logPrior+=log(*(REAL8 *)LALInferenceGetVariable(params,"flow"));
1291  }
1292 
1293  if(LALInferenceCheckVariable(params,"distance"))
1294  logPrior-=log(*(REAL8 *)LALInferenceGetVariable(params,"distance"));
1295  if(LALInferenceCheckVariable(params,"theta_jn"))
1296  logPrior+=log(fabs(sin(*(REAL8 *)LALInferenceGetVariable(params,"theta_jn"))));
1297  if(LALInferenceCheckVariable(params,"declination"))
1298  {
1299  /* Check that this is not an output variable */
1300  if(LALInferenceGetVariableVaryType(params,"declination")==LALINFERENCE_PARAM_LINEAR)
1301  logPrior+=log(fabs(cos(*(REAL8 *)LALInferenceGetVariable(params,"declination"))));
1302  }
1303  if(LALInferenceCheckVariable(params,"tilt_spin1"))
1304  {
1305  LALInferenceParamVaryType vtype=LALInferenceGetVariableVaryType(params,"tilt_spin1");
1306  if(vtype!=LALINFERENCE_PARAM_FIXED && vtype!=LALINFERENCE_PARAM_OUTPUT)
1307  {
1308  logPrior+=log(fabs(sin(*(REAL8 *)LALInferenceGetVariable(params,"tilt_spin1"))));
1309  }
1310  }
1311  if(LALInferenceCheckVariable(params,"tilt_spin2"))
1312  {
1313  LALInferenceParamVaryType vtype=LALInferenceGetVariableVaryType(params,"tilt_spin2");
1314  if(vtype!=LALINFERENCE_PARAM_FIXED && vtype!=LALINFERENCE_PARAM_OUTPUT)
1315  {
1316  logPrior+=log(fabs(sin(*(REAL8 *)LALInferenceGetVariable(params,"tilt_spin2"))));
1317  }
1318  }
1319  /*priors uniform in the individual masses. Not taking into account if mtot_max < m1_max+m2_max */
1320  if(LALInferenceCheckVariable(params,"eta")||LALInferenceCheckVariable(params,"q")) {
1321  if(LALInferenceCheckVariable(params,"logmc")) {
1322  logmc=*(REAL8 *)LALInferenceGetVariable(params,"logmc");
1323  if(LALInferenceCheckVariable(params,"q")) {
1324  q=*(REAL8 *)LALInferenceGetVariable(params,"q");
1325  LALInferenceMcQ2Masses(exp(logmc),q,&m1,&m2);
1326  logPrior+=log(m1*m1);
1327  } else {
1328  eta=*(REAL8 *)LALInferenceGetVariable(params,"eta");
1329  LALInferenceMcEta2Masses(exp(logmc),eta,&m1,&m2);
1330  logPrior+=log(((m1+m2)*(m1+m2)*(m1+m2))/(m1-m2));
1331  }
1332  /*careful using LALInferenceMcEta2Masses, it returns m1>=m2*/
1333  } else if(LALInferenceCheckVariable(params,"chirpmass")) {
1334  mc=*(REAL8 *)LALInferenceGetVariable(params,"chirpmass");
1335  if(LALInferenceCheckVariable(params,"q")) {
1336  q=*(REAL8 *)LALInferenceGetVariable(params,"q");
1337  LALInferenceMcQ2Masses(mc,q,&m1,&m2);
1338  logPrior+=log(m1*m1/mc);
1339  } else {
1340  eta=*(REAL8 *)LALInferenceGetVariable(params,"eta");
1341  LALInferenceMcEta2Masses(mc,eta,&m1,&m2);
1342  logPrior+=log(((m1+m2)*(m1+m2))/((m1-m2)*pow(eta,3.0/5.0)));
1343  }
1344  }
1345  }
1346 
1347  if(LALInferenceCheckVariable(priorParams,"MTotMax"))
1348  if(*(REAL8 *)LALInferenceGetVariable(priorParams,"MTotMax") < m1+m2)
1349  return -INFINITY;
1350 
1351  if(LALInferenceCheckVariable(priorParams,"MTotMin"))
1352  if(*(REAL8 *)LALInferenceGetVariable(priorParams,"MTotMin") > m1+m2)
1353  return -INFINITY;
1354  /* Check for individual mass priors */
1355  if(LALInferenceCheckVariable(priorParams,"mass1_min"))
1356  if(LALInferenceGetREAL8Variable(priorParams,"mass1_min") > m1)
1357  return -INFINITY;
1358  if(LALInferenceCheckVariable(priorParams,"mass1_max"))
1359  if(LALInferenceGetREAL8Variable(priorParams,"mass1_max") < m1)
1360  return -INFINITY;
1361  if(LALInferenceCheckVariable(priorParams,"mass2_min"))
1362  if(LALInferenceGetREAL8Variable(priorParams,"mass2_min") > m2)
1363  return -INFINITY;
1364  if(LALInferenceCheckVariable(priorParams,"mass2_max"))
1365  if(LALInferenceGetREAL8Variable(priorParams,"mass2_max") < m2)
1366  return -INFINITY;
1367 
1368  //PSD priors are Gaussian
1369  if(LALInferenceCheckVariable(params, "psdscale"))
1370  {
1371  UINT4 i;
1372  UINT4 j;
1373 
1374  REAL8 var;
1375  REAL8 mean = 1.0;
1376  REAL8 prior= 0.0;
1377  UINT4 psdGaussianPrior;
1378 
1379  REAL8Vector *sigma = *((REAL8Vector **)LALInferenceGetVariable(priorParams, "psdsigma"));
1380  gsl_matrix *nparams = *((gsl_matrix **)LALInferenceGetVariable(params,"psdscale"));
1381 
1382  min=*(REAL8 *)LALInferenceGetVariable(priorParams,"psdrange_min");
1383  max=*(REAL8 *)LALInferenceGetVariable(priorParams,"psdrange_max");
1384 
1385  psdGaussianPrior = *(UINT4 *)LALInferenceGetVariable(priorParams,"psdGaussianPrior");
1386 
1387  for(i=0; i<(UINT4)nparams->size1; i++)
1388  {
1389  for(j=0; j<(UINT4)nparams->size2; j++)
1390  {
1391  var = sigma->data[j]*sigma->data[j];
1392  val = gsl_matrix_get(nparams,i,j);
1393  //reject prior
1394  if(val < min || val > max) return -INFINITY;
1395  else if(psdGaussianPrior)prior += -0.5*( (mean-val)*(mean-val)/var + log(2.0*LAL_PI*var) );
1396  }
1397  }
1398  logPrior+=prior;
1399  }
1400 
1401  /* Calibration parameters */
1402  //logPrior += LALInferenceSplineCalibrationPrior(runState, params);
1403  logPrior += LALInferenceConstantCalibrationPrior(runState, params);
1404  return(logPrior);
1405 }
1406 
1407 /* Convert the hypercube parameter to physical parameters, for the non-spinning inspiral signal case */
1408 UINT4 LALInferenceInspiralSkyLocCubeToPrior(LALInferenceRunState *runState, LALInferenceVariables *params, UNUSED LALInferenceModel *model, double *Cube, void *context)
1409 {
1410  LALInferenceVariableItem *item;
1411  LALInferenceVariables *priorParams=runState->priorArgs;
1412 
1413  char **info = (char **)context;
1414  char *timeID = &info[2][0];
1415  REAL8 min=-INFINITY, max=INFINITY;
1416  int i = 0;
1417 
1418  INT4 SKY_FRAME=0;
1419  if(LALInferenceCheckVariable(params,"SKY_FRAME"))
1420  SKY_FRAME=*(INT4 *)LALInferenceGetVariable(params,"SKY_FRAME");
1421  double azimuth, cosalpha, lat, longitude, t0, tc;
1422 
1423  if(SKY_FRAME==1)
1424  {
1425  // detector frame azimuth
1426  if (LALInferenceCheckVariable(params, "azimuth"))
1427  {
1428  item = LALInferenceGetItem(params, "azimuth");
1429  if (item->vary != LALINFERENCE_PARAM_FIXED)
1430  {
1431  LALInferenceGetMinMaxPrior(runState->priorArgs, "azimuth", (void *)&min, (void *)&max);
1432  azimuth = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1433  LALInferenceSetVariable(params, "azimuth", &azimuth);
1434  i++;
1435  }
1436  }
1437 
1438  // detector frame cosalpha
1439  if (LALInferenceCheckVariable(params, "cosalpha"))
1440  {
1441  item = LALInferenceGetItem(params, "cosalpha");
1442  if (item->vary != LALINFERENCE_PARAM_FIXED)
1443  {
1444  LALInferenceGetMinMaxPrior(runState->priorArgs, "cosalpha", (void *)&min, (void *)&max);
1445  cosalpha = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1446  LALInferenceSetVariable(params, "cosalpha", &cosalpha);
1447  i++;
1448  }
1449  }
1450  }
1451  else
1452  {
1453  // latitude
1454  if (LALInferenceCheckVariable(params, "declination"))
1455  {
1456  item = LALInferenceGetItem(params, "declination");
1457  if(item->vary != LALINFERENCE_PARAM_FIXED)
1458  {
1459  lat = asin(2.0 * Cube[i] - 1.0);
1460  LALInferenceSetVariable(params, "declination", &lat);
1461  i++;
1462  }
1463  }
1464 
1465  // longitude
1466  if (LALInferenceCheckVariable(params, "rightascension"))
1467  {
1468  item = LALInferenceGetItem(params, "rightascension");
1469  if(item->vary != LALINFERENCE_PARAM_FIXED)
1470  {
1471  LALInferenceGetMinMaxPrior(runState->priorArgs, "rightascension", (void *)&min, (void *)&max);
1472  longitude = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1473  LALInferenceSetVariable(params, "rightascension", &longitude);
1474  i++;
1475  }
1476  }
1477  }
1478 
1479  // phi
1480  if (LALInferenceCheckVariable(params, "phase"))
1481  {
1482  item = LALInferenceGetItem(params, "phase");
1483  if(item->vary != LALINFERENCE_PARAM_FIXED)
1484  {
1485  LALInferenceGetMinMaxPrior(runState->priorArgs, "phase", (void *)&min, (void *)&max);
1486  double phi = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1487  LALInferenceSetVariable(params, "phase", &phi);
1488  i++;
1489  }
1490  }
1491 
1492  // psi
1493  if (LALInferenceCheckVariable(params, "polarisation"))
1494  {
1495  item = LALInferenceGetItem(params, "polarisation");
1496  if(item->vary != LALINFERENCE_PARAM_FIXED)
1497  {
1498  LALInferenceGetMinMaxPrior(runState->priorArgs, "polarisation", (void *)&min, (void *)&max);
1499  double psi = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1500  LALInferenceSetVariable(params, "polarisation", &psi);
1501  i++;
1502  }
1503  }
1504 
1505  if(SKY_FRAME==1)
1506  {
1507  // detector frame t0
1508  if (LALInferenceCheckVariable(params, "t0"))
1509  {
1510  item = LALInferenceGetItem(params, "t0");
1511  if(item->vary != LALINFERENCE_PARAM_FIXED)
1512  {
1513  LALInferenceGetMinMaxPrior(runState->priorArgs, "t0", (void *)&min, (void *)&max);
1514  t0 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1515  LALInferenceSetVariable(params, "t0", &t0);
1516  sprintf(timeID,"%d",i);
1517  i++;
1518  }
1519  }
1520  }
1521  else
1522  {
1523  // time
1524  if (LALInferenceCheckVariable(params, "time"))
1525  {
1526  item = LALInferenceGetItem(params, "time");
1527  if(item->vary != LALINFERENCE_PARAM_FIXED)
1528  {
1529  LALInferenceGetMinMaxPrior(runState->priorArgs, "time", (void *)&min, (void *)&max);
1530  tc = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1531  LALInferenceSetVariable(params, "time", &tc);
1532  sprintf(timeID,"%d",i);
1533  i++;
1534  }
1535  }
1536  }
1537 
1538  // check if mchirp is fixed
1539  double mc = 0.0, eta = 0.0, q = 0.0, m1 = 0.0, m2 = 0.0, m = 0.0;
1540  if( LALInferenceCheckVariable(params,"logmc") )
1541  {
1542  item = LALInferenceGetItem(params, "logmc");
1543  if(item->vary == LALINFERENCE_PARAM_FIXED) mc = exp(*(REAL8 *)LALInferenceGetVariable(params, "logmc"));
1544  }
1545  else if( LALInferenceCheckVariable(params,"chirpmass") )
1546  {
1547  item = LALInferenceGetItem(params, "chirpmass");
1548  if(item->vary == LALINFERENCE_PARAM_FIXED) mc = *(REAL8 *)LALInferenceGetVariable(params, "chirpmass");
1549  }
1550 
1551  // check if eta is fixed
1552  if( LALInferenceCheckVariable(params,"eta") )
1553  {
1554  item = LALInferenceGetItem(params, "eta");
1555  if(item->vary == LALINFERENCE_PARAM_FIXED)
1556  {
1557  eta = *(REAL8 *)LALInferenceGetVariable(params, "eta");
1558  if( mc != 0.0 ) LALInferenceMcEta2Masses(mc,eta,&m1,&m2);
1559  }
1560  }
1561  else if( LALInferenceCheckVariable(params,"q") )
1562  {
1563  item = LALInferenceGetItem(params, "q");
1564  if(item->vary == LALINFERENCE_PARAM_FIXED)
1565  {
1566  q = *(REAL8 *)LALInferenceGetVariable(params, "q");
1567  if( mc != 0.0 ) LALInferenceMcQ2Masses(mc,q,&m1,&m2);
1568  }
1569  }
1570 
1571  //m1 & m2
1572  if( m1 == 0.0 && m2 == 0.0 )
1573  {
1574  REAL8 m1_min = *(REAL8 *)LALInferenceGetVariable(priorParams,"mass1_min");
1575  REAL8 m1_max = *(REAL8 *)LALInferenceGetVariable(priorParams,"mass1_max");
1576  REAL8 m2_min = *(REAL8 *)LALInferenceGetVariable(priorParams,"mass2_min");
1577  REAL8 m2_max = *(REAL8 *)LALInferenceGetVariable(priorParams,"mass2_max");
1578  if( m1_min == m1_max && m2_min == m2_max)
1579  {
1580  m1 = m1_min;
1581  m2 = m2_min;
1582  m = m1 + m2;
1583  eta = m1 * m2 / (m*m);
1584  mc = pow(eta,0.6) * m;
1585  q = m2 / m1; // asymmetric mass ratio, m1 >= m2
1586  }
1587  else
1588  {
1589  m1 = LALInferenceCubeToFlatPrior(Cube[i], m1_min, m1_max);
1590  m2 = LALInferenceCubeToFlatPrior(Cube[i+1], m2_min, m2_max);
1591  if(m1<m2)
1592  {
1593  double temp = m2;
1594  m2 = m1;
1595  m1 = temp;
1596  }
1597 
1598  m = m1 + m2;
1599  eta = m1 * m2 / (m*m);
1600  mc = pow(eta,0.6) * m;
1601  q = m2 / m1; // asymmetric mass ratio, m1 >= m2
1602  i++;
1603  i++;
1604  }
1605 
1606  // chirp mass and eta/q
1607  if(LALInferenceCheckVariable(params,"eta")||LALInferenceCheckVariable(params,"q"))
1608  {
1609  if(LALInferenceCheckVariable(params,"logmc"))
1610  {
1611  double logmc = log(mc);
1612  LALInferenceSetVariable(params, "logmc", &logmc);
1613  }
1614  else if(LALInferenceCheckVariable(params,"chirpmass"))
1615  {
1616  LALInferenceSetVariable(params, "chirpmass", &mc);
1617  }
1618 
1619  if(LALInferenceCheckVariable(params,"q"))
1620  LALInferenceSetVariable(params, "q", &q);
1621  else if(LALInferenceCheckVariable(params,"eta"))
1622  LALInferenceSetVariable(params, "eta", &eta);
1623  }
1624  }
1625 
1626  // distance
1627  double dist;
1628  if( LALInferenceCheckVariable(params,"logdistance") )
1629  {
1630  item = LALInferenceGetItem(params, "logdistance");
1631  if(item->vary != LALINFERENCE_PARAM_FIXED)
1632  {
1633  LALInferenceGetMinMaxPrior(runState->priorArgs, "logdistance", (void *)&min, (void *)&max);
1634  min = exp(min); max = exp(max);
1635  dist = LALInferenceCubeToLogFlatPrior(Cube[i], min, max);
1636  double logdist = log(dist);
1637  LALInferenceSetVariable(params, "logdistance", &logdist);
1638  i++;
1639  }
1640  }
1641  else if( LALInferenceCheckVariable(params,"distance") )
1642  {
1643  item = LALInferenceGetItem(params, "distance");
1644  if(item->vary != LALINFERENCE_PARAM_FIXED)
1645  {
1646  LALInferenceGetMinMaxPrior(runState->priorArgs, "distance", (void *)&min, (void *)&max);
1647  dist = LALInferenceCubeToLogFlatPrior(Cube[i], min, max);
1648  LALInferenceSetVariable(params, "distance", &dist);
1649  i++;
1650  }
1651  }
1652 
1653  // a_spin1
1654  if(LALInferenceCheckVariable(params,"a_spin1"))
1655  {
1656  item = LALInferenceGetItem(params, "a_spin1");
1657  if(item->vary != LALINFERENCE_PARAM_FIXED)
1658  {
1659  LALInferenceGetMinMaxPrior(runState->priorArgs, "a_spin1", (void *)&min, (void *)&max);
1660  double a_spin1 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1661  LALInferenceSetVariable(params, "a_spin1", &a_spin1);
1662  i++;
1663  }
1664  }
1665 
1666  // a_spin2
1667  if(LALInferenceCheckVariable(params,"a_spin2"))
1668  {
1669  item = LALInferenceGetItem(params, "a_spin2");
1670  if(item->vary != LALINFERENCE_PARAM_FIXED)
1671  {
1672  LALInferenceGetMinMaxPrior(runState->priorArgs, "a_spin2", (void *)&min, (void *)&max);
1673  double a_spin2 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1674  LALInferenceSetVariable(params, "a_spin2", &a_spin2);
1675  i++;
1676  }
1677  }
1678 
1679  // theta_JN for system-frame parameters
1680  if(LALInferenceCheckVariable(params,"costheta_jn"))
1681  {
1682  item = LALInferenceGetItem(params, "costheta_jn");
1683  if(item->vary != LALINFERENCE_PARAM_FIXED)
1684  {
1685  LALInferenceGetMinMaxPrior(runState->priorArgs, "costheta_jn", (void *)&min, (void *)&max);
1686  double costheta_JN = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1687  LALInferenceSetVariable(params, "costheta_jn", &costheta_JN);
1688  i++;
1689  }
1690  }
1691 
1692  // phi_JL for system-frame parameters
1693  if(LALInferenceCheckVariable(params,"phi_jl"))
1694  {
1695  item = LALInferenceGetItem(params, "phi_jl");
1696  if(item->vary != LALINFERENCE_PARAM_FIXED)
1697  {
1698  LALInferenceGetMinMaxPrior(runState->priorArgs, "phi_jl", (void *)&min, (void *)&max);
1699  double phi_JL = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1700  LALInferenceSetVariable(params, "phi_jl", &phi_JL);
1701  i++;
1702  }
1703  }
1704 
1705  // tilt of spin 1 for system-frame parameters
1706  if(LALInferenceCheckVariable(params,"tilt_spin1"))
1707  {
1708  item = LALInferenceGetItem(params, "tilt_spin1");
1709  if(item->vary != LALINFERENCE_PARAM_FIXED)
1710  {
1711  LALInferenceGetMinMaxPrior(runState->priorArgs, "tilt_spin1", (void *)&min, (void *)&max);
1712  double tilt_spin1 = LALInferenceCubeToSinPrior(Cube[i], min, max);
1713  LALInferenceSetVariable(params, "tilt_spin1", &tilt_spin1);
1714  i++;
1715  }
1716  }
1717 
1718  // tilt of spin 2 for system-frame parameters
1719  if(LALInferenceCheckVariable(params,"tilt_spin2"))
1720  {
1721  item = LALInferenceGetItem(params, "tilt_spin2");
1722  if(item->vary != LALINFERENCE_PARAM_FIXED)
1723  {
1724  LALInferenceGetMinMaxPrior(runState->priorArgs, "tilt_spin2", (void *)&min, (void *)&max);
1725  double tilt_spin2 = LALInferenceCubeToSinPrior(Cube[i], min, max);
1726  LALInferenceSetVariable(params, "tilt_spin2", &tilt_spin2);
1727  i++;
1728  }
1729  }
1730 
1731  // phi12 for system-frame parameters
1732  if(LALInferenceCheckVariable(params,"phi12"))
1733  {
1734  item = LALInferenceGetItem(params, "phi12");
1735  if(item->vary != LALINFERENCE_PARAM_FIXED)
1736  {
1737  LALInferenceGetMinMaxPrior(runState->priorArgs, "phi12", (void *)&min, (void *)&max);
1738  double phi12 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
1739  LALInferenceSetVariable(params, "phi12", &phi12);
1740  i++;
1741  }
1742  }
1743 
1744  INT4 ScaleTest = LALInferenceCubeToPSDScaleParams(priorParams, params, &i, Cube, context);
1745  UINT4 ConstCalib = LALInferenceCubeToConstantCalibrationPrior(runState, params, &i, Cube, context);
1746  //UINT4 SplineCalib = LALInferenceCubeToSplineCalibrationPrior(runState, params, &i, Cube, context);
1747 
1748  /* Check boundaries */
1749  if (ScaleTest==0 || ConstCalib==0 /* || SplineCalib==0 */) return 0;
1750  item=params->head;
1751  for(;item;item=item->next)
1752  {
1753  if(item->vary==LALINFERENCE_PARAM_FIXED || item->vary==LALINFERENCE_PARAM_OUTPUT)
1754  continue;
1755  else
1756  {
1757  LALInferenceGetMinMaxPrior(priorParams, item->name, (void *)&min, (void *)&max);
1758  if(*(REAL8 *) item->value < min || *(REAL8 *)item->value > max) return 0 ;
1759  }
1760  }
1761 
1762  /* Check for real mass values */
1763  if(LALInferenceCheckVariable(params,"logmc"))
1764  if(isnan(*(REAL8 *)LALInferenceGetVariable(params,"logmc")))
1765  return 0;
1766  if(LALInferenceCheckVariable(params,"chirpmass"))
1767  if(isnan(*(REAL8 *)LALInferenceGetVariable(params,"chirpmass")))
1768  return 0;
1769  if(LALInferenceCheckVariable(params,"eta"))
1770  if(isnan(*(REAL8 *)LALInferenceGetVariable(params,"eta"))
1771  ||*(REAL8 *)LALInferenceGetVariable(params,"eta") < 0.0
1772  || *(REAL8 *)LALInferenceGetVariable(params,"eta") > 0.25)
1773  return 0;
1774  if(LALInferenceCheckVariable(params,"q"))
1775  if(isnan(*(REAL8 *)LALInferenceGetVariable(params,"q"))
1776  ||*(REAL8 *)LALInferenceGetVariable(params,"q") < 0.0
1777  || *(REAL8 *)LALInferenceGetVariable(params,"q") > 1.0)
1778  return 0;
1779 
1780  if(LALInferenceCheckVariable(priorParams,"MTotMin"))
1781  if(*(REAL8 *)LALInferenceGetVariable(priorParams,"MTotMin") > m1+m2)
1782  return 0;
1783 
1784  if(LALInferenceCheckVariable(priorParams,"MTotMax"))
1785  if(*(REAL8 *)LALInferenceGetVariable(priorParams,"MTotMax") < m1+m2)
1786  return 0;
1787  /* Check for individual mass priors */
1788  if(LALInferenceCheckVariable(priorParams,"mass1_min"))
1789  if(LALInferenceGetREAL8Variable(priorParams,"mass1_min") > m1)
1790  return 0;
1791  if(LALInferenceCheckVariable(priorParams,"mass1_max"))
1792  if(LALInferenceGetREAL8Variable(priorParams,"mass1_max") < m1)
1793  return 0;
1794  if(LALInferenceCheckVariable(priorParams,"mass2_min"))
1795  if(LALInferenceGetREAL8Variable(priorParams,"mass2_min") > m2)
1796  return 0;
1797  if(LALInferenceCheckVariable(priorParams,"mass2_max"))
1798  if(LALInferenceGetREAL8Variable(priorParams,"mass2_max") < m2)
1799  return 0;
1800 
1801  return 1;
1802 }
1803 
1804 
1805 typedef struct {
1806  double M1;
1807  double McMin;
1808  double McMax;
1809  double massRatioMin;
1810  double massRatioMax;
1811 } innerData;
1812 
1813 static double qInnerIntegrand(double M2, void *viData) {
1814  innerData *iData = (innerData *)viData;
1815  double Mc = pow(M2*iData->M1, 3.0/5.0)/pow(M2+iData->M1, 1.0/5.0);
1816  double q = M2/iData->M1;
1817  if (Mc < iData->McMin || Mc > iData->McMax || q < iData->massRatioMin || q > iData->massRatioMax) {
1818  return 0.0;
1819  } else {
1820  return pow(Mc, -11.0/6.0);
1821  }
1822 }
1823 
1824 #define LALINFERENCE_PRIOR_SQR(x) ((x)*(x))
1825 
1826 static double etaInnerIntegrand(double M2, void *viData) {
1827  innerData *iData = (innerData *)viData;
1828  double Mc = pow(M2*iData->M1, 3.0/5.0)/pow(M2+iData->M1, 1.0/5.0);
1829  double eta = M2*iData->M1/LALINFERENCE_PRIOR_SQR(M2+iData->M1);
1830  if (Mc < iData->McMin || Mc > iData->McMax || eta < iData->massRatioMin || eta > iData->massRatioMax) {
1831  return 0.0;
1832  } else {
1833  return pow(Mc, -11.0/6.0);
1834  }
1835 }
1836 
1837 #undef LALINFERENCE_PRIOR_SQR
1838 
1839 typedef struct {
1840  gsl_integration_workspace *wsInner;
1841  size_t wsInnerSize;
1842  double McMin;
1843  double McMax;
1844  double massRatioMin;
1845  double massRatioMax;
1846  double MTotMax;
1847  double MMin;
1848  double epsabs;
1849  double epsrel;
1850  gsl_function innerIntegrand;
1851 } outerData;
1852 
1853 #define LALINFERENCE_PRIOR_MIN(x, y) ((x) < (y) ? (x) : (y))
1854 
1855 static double outerIntegrand(double M1, void *voData) {
1856 
1857  outerData *oData = (outerData *)voData;
1858  gsl_function f;
1859  innerData iData;
1860  double result, err;
1861 
1862  iData.M1 = M1;
1863  iData.McMin = oData->McMin;
1864  iData.McMax = oData->McMax;
1865  iData.massRatioMin = oData->massRatioMin;
1866  iData.massRatioMax = oData->massRatioMax;
1867 
1868  f.function=(oData->innerIntegrand.function);
1869 
1870  f.params = &iData;
1871 
1872  gsl_integration_qag(&f, oData->MMin, LALINFERENCE_PRIOR_MIN(M1, oData->MTotMax-M1), oData->epsabs, oData->epsrel,
1873  oData->wsInnerSize, GSL_INTEG_GAUSS61, oData->wsInner, &result, &err);
1874 
1875  return result;
1876 }
1877 
1878 #undef LALINFERENCE_PRIOR_MIN
1879 
1880 REAL8 LALInferenceComputePriorMassNorm(const double MMin, const double MMax, const double MTotMax,
1881  const double McMin, const double McMax,
1882  const double massRatioMin, const double massRatioMax, const char *massRatioName) {
1883  const double epsabs = 1e-8;
1884  const double epsrel = 1e-8;
1885  const size_t wsSize = 10000;
1886  double result, err;
1887  outerData oData;
1888  gsl_function f;
1889 
1890  if(!massRatioName)
1891  XLAL_ERROR_REAL8(XLAL_EFAULT, "Null arguments received.");
1892  else if(!strcmp(massRatioName,"q"))
1893  oData.innerIntegrand.function = &qInnerIntegrand;
1894  else if(!strcmp(massRatioName,"eta"))
1895  oData.innerIntegrand.function = &etaInnerIntegrand;
1896  else
1897  XLAL_ERROR_REAL8(XLAL_ENAME, "Invalid mass ratio name specified");
1898 
1899  // Disable GSL error reporting in favour of XLAL (the integration routines are liable to fail).
1900  gsl_error_handler_t *oldHandler = gsl_set_error_handler_off();
1901 
1902  gsl_integration_workspace *wsOuter = gsl_integration_workspace_alloc(wsSize);
1903  gsl_integration_workspace *wsInner = gsl_integration_workspace_alloc(wsSize);
1904 
1905  oData.wsInnerSize = wsSize;
1906  oData.wsInner = wsInner;
1907  oData.McMin = McMin;
1908  oData.McMax = McMax;
1909  oData.massRatioMin = massRatioMin;
1910  oData.massRatioMax = massRatioMax;
1911  oData.MTotMax = MTotMax;
1912  oData.epsabs = epsabs;
1913  oData.epsrel = epsrel;
1914  oData.MMin = MMin;
1915 
1916  f.function = &outerIntegrand;
1917  f.params = &oData;
1918 
1919  int status = gsl_integration_qag(&f, MMin, MMax, epsabs, epsrel, wsSize, GSL_INTEG_GAUSS61, wsOuter,
1920  &result, &err);
1921 
1922  if (status)
1923  XLAL_ERROR_REAL8(XLAL_EFUNC | XLAL_EDATA, "Bad data; GSL integration failed.");
1924 
1925  gsl_set_error_handler(oldHandler);
1926 
1927  gsl_integration_workspace_free(wsOuter);
1928  gsl_integration_workspace_free(wsInner);
1929 
1930  return result;
1931 }
1932 
1933 /* Function to add the min and max values for the prior onto the priorArgs */
1934 void LALInferenceAddMinMaxPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *min, REAL8 *max, LALInferenceVariableType type){
1935  if (*min >= *max)
1936  XLAL_ERROR_VOID(XLAL_EINVAL, "Minimum must be less than maximum, but %f >= %f.", *min, *max);
1937 
1938  char minName[VARNAME_MAX];
1939  char maxName[VARNAME_MAX];
1940 
1941  sprintf(minName,"%s_min",name);
1942  sprintf(maxName,"%s_max",name);
1943 
1944  LALInferenceAddVariable(priorArgs,minName,min,type,LALINFERENCE_PARAM_FIXED);
1945  LALInferenceAddVariable(priorArgs,maxName,max,type,LALINFERENCE_PARAM_FIXED);
1946  return;
1947 }
1948 
1949 /* Function to remove the min and max values for the prior onto the priorArgs */
1950 void LALInferenceRemoveMinMaxPrior(LALInferenceVariables *priorArgs, const char *name){
1951  char minName[VARNAME_MAX];
1952  char maxName[VARNAME_MAX];
1953 
1954  sprintf(minName,"%s_min",name);
1955  sprintf(maxName,"%s_max",name);
1956 
1957  LALInferenceRemoveVariable(priorArgs, minName);
1958  LALInferenceRemoveVariable(priorArgs, maxName);
1959  return;
1960 }
1961 
1962 /* Check for a min/max prior set */
1963 int LALInferenceCheckMinMaxPrior(LALInferenceVariables *priorArgs, const char *name)
1964 {
1965  char minName[VARNAME_MAX+4];
1966  char maxName[VARNAME_MAX+4];
1967  sprintf(minName,"%s_min",name);
1968  sprintf(maxName,"%s_max",name);
1969 
1970  return (LALInferenceCheckVariable(priorArgs,minName) && LALInferenceCheckVariable(priorArgs,maxName));
1971 }
1972 
1973 /* Get the min and max values of the prior from the priorArgs list, given a name */
1974 void LALInferenceGetMinMaxPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *min, REAL8 *max)
1975 {
1976  char minName[VARNAME_MAX+4];
1977  char maxName[VARNAME_MAX+4];
1978  void *ptr=NULL;
1979  sprintf(minName,"%s_min",name);
1980  sprintf(maxName,"%s_max",name);
1981 
1982  ptr=LALInferenceGetVariable(priorArgs,minName);
1983  if(ptr) *min=*(REAL8*)ptr;
1984  else XLAL_ERROR_VOID(XLAL_EFAILED);
1985  ptr=LALInferenceGetVariable(priorArgs,maxName);
1986  if(ptr) *max=*(REAL8*)ptr;
1987  else XLAL_ERROR_VOID(XLAL_EFAILED);
1988  return;
1989 }
1990 
1991 
1992 /* Check for a Gaussian Prior of the standard form */
1993 int LALInferenceCheckGaussianPrior(LALInferenceVariables *priorArgs, const char *name)
1994 {
1995  char meanName[VARNAME_MAX+14];
1996  char sigmaName[VARNAME_MAX+15];
1997  sprintf(meanName,"%s_gaussian_mean",name);
1998  sprintf(sigmaName,"%s_gaussian_sigma",name);
1999  return (LALInferenceCheckVariable(priorArgs,meanName) && LALInferenceCheckVariable(priorArgs,sigmaName));
2000 }
2001 
2002 /* Function to add the min and max values for the prior onto the priorArgs */
2003 void LALInferenceAddGaussianPrior( LALInferenceVariables *priorArgs,
2004  const char *name, REAL8 *mu, REAL8 *sigma,
2005  LALInferenceVariableType type ){
2006  char meanName[VARNAME_MAX+14];
2007  char sigmaName[VARNAME_MAX+15];
2008 
2009  sprintf(meanName,"%s_gaussian_mean",name);
2010  sprintf(sigmaName,"%s_gaussian_sigma",name);
2011 
2012  LALInferenceAddVariable(priorArgs,meanName,mu,type,LALINFERENCE_PARAM_FIXED);
2013  LALInferenceAddVariable(priorArgs,sigmaName,sigma,type,LALINFERENCE_PARAM_FIXED);
2014  return;
2015 }
2016 
2017 /* Function to remove the min and max values for the prior onto the priorArgs */
2018 void LALInferenceRemoveGaussianPrior(LALInferenceVariables *priorArgs, const char *name){
2019  char meanName[VARNAME_MAX];
2020  char sigmaName[VARNAME_MAX];
2021 
2022  sprintf(meanName,"%s_gaussian_mean",name);
2023  sprintf(sigmaName,"%s_gaussian_sigma",name);
2024 
2025  LALInferenceRemoveVariable(priorArgs, meanName);
2026  LALInferenceRemoveVariable(priorArgs, sigmaName);
2027  return;
2028 }
2029 
2030 /* Get the min and max values of the prior from the priorArgs list, given a name
2031 */
2032 void LALInferenceGetGaussianPrior(LALInferenceVariables *priorArgs,
2033  const char *name, REAL8 *mu, REAL8 *sigma)
2034 {
2035  char meanName[VARNAME_MAX];
2036  char sigmaName[VARNAME_MAX];
2037  void *ptr=NULL;
2038 
2039  sprintf(meanName,"%s_gaussian_mean",name);
2040  sprintf(sigmaName,"%s_gaussian_sigma",name);
2041 
2042  ptr = LALInferenceGetVariable(priorArgs, meanName);
2043  if ( ptr ) *mu = *(REAL8*)ptr;
2044  else XLAL_ERROR_VOID(XLAL_EFAILED);
2045 
2046  ptr = LALInferenceGetVariable(priorArgs, sigmaName);
2047  if ( ptr ) *sigma = *(REAL8*)ptr;
2048  else XLAL_ERROR_VOID(XLAL_EFAILED);
2049 
2050  return;
2051 }
2052 
2053 /* Function to add a correlation matrix to the prior onto the priorArgs */
2054 void LALInferenceAddCorrelatedPrior(LALInferenceVariables *priorArgs,
2055  const char *name, gsl_matrix **cor,
2056  REAL8 *mu, REAL8 *sigma, UINT4 *idx){
2057  char corName[VARNAME_MAX];
2058  char invName[VARNAME_MAX];
2059  char meanName[VARNAME_MAX];
2060  char sigmaName[VARNAME_MAX];
2061  char idxName[VARNAME_MAX];
2062 
2063  sprintf(corName,"correlation_matrix");
2064  sprintf(invName,"inverse_correlation_matrix");
2065 
2066  /* add correlation matrix if not already added */
2067  if ( !LALInferenceCheckVariable( priorArgs, corName ) ){
2068  LALInferenceAddVariable(priorArgs, corName, cor, LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_FIXED);
2069 
2070  /* get the matrix we've just added */
2071  gsl_matrix *thiscor = NULL, *usecor = NULL;
2072  thiscor = *(gsl_matrix **)LALInferenceGetVariable( priorArgs, corName );
2073  usecor = gsl_matrix_alloc( thiscor->size1, thiscor->size2 );
2074  XLAL_CALLGSL( gsl_matrix_memcpy( usecor, thiscor ) );
2075 
2076  gsl_matrix *invcor = gsl_matrix_alloc( usecor->size1, usecor->size2 );
2077  gsl_permutation *p = gsl_permutation_alloc( usecor->size1 );
2078  INT4 s;
2079 
2080  /* check correlation matrix is positive definite */
2081  if( !LALInferenceCheckPositiveDefinite( usecor, usecor->size1 ) ){
2082  XLAL_ERROR_VOID( XLAL_EFUNC | XLAL_EINVAL, "Error... matrix is not positive definite!" );
2083  }
2084 
2085  /* invert correlation matrix */
2086  XLAL_CALLGSL( gsl_linalg_LU_decomp( usecor, p, &s ) );
2087  XLAL_CALLGSL( gsl_linalg_LU_invert( usecor, p, invcor ) );
2088  XLAL_CALLGSL( gsl_permutation_free( p ) );
2089  XLAL_CALLGSL( gsl_matrix_free( usecor) );
2090 
2091  LALInferenceAddVariable(priorArgs, invName, &invcor, LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_FIXED);
2092  }
2093 
2094  sprintf(meanName, "%s_correlation_mean", name);
2095  sprintf(sigmaName, "%s_correlation_sigma", name);
2096  sprintf(idxName,"%s_index",name);
2097 
2098  LALInferenceAddVariable(priorArgs, meanName, mu, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
2099  LALInferenceAddVariable(priorArgs, sigmaName, sigma, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
2100  LALInferenceAddVariable(priorArgs, idxName, idx, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);
2101 
2102  return;
2103 }
2104 
2105 /* Get the correlation matrix and parameter index */
2106 void LALInferenceGetCorrelatedPrior(LALInferenceVariables *priorArgs,
2107  const char *name, gsl_matrix **cor, gsl_matrix **invcor,
2108  REAL8 *mu, REAL8 *sigma, UINT4 *idx){
2109  char corName[VARNAME_MAX];
2110  char invName[VARNAME_MAX];
2111  char meanName[VARNAME_MAX];
2112  char sigmaName[VARNAME_MAX];
2113  char idxName[VARNAME_MAX];
2114  void *ptr = NULL;
2115 
2116  sprintf(corName, "correlation_matrix");
2117  sprintf(invName,"inverse_correlation_matrix");
2118  sprintf(meanName, "%s_correlation_mean", name);
2119  sprintf(sigmaName, "%s_correlation_sigma", name);
2120  sprintf(idxName,"%s_index",name);
2121 
2122  ptr = LALInferenceGetVariable(priorArgs, corName);
2123  if ( ptr ) *cor = *(gsl_matrix **)ptr;
2124  else XLAL_ERROR_VOID(XLAL_EFAILED);
2125 
2126  ptr = LALInferenceGetVariable(priorArgs, invName);
2127  if ( ptr ) *invcor = *(gsl_matrix **)ptr;
2128  else XLAL_ERROR_VOID(XLAL_EFAILED);
2129 
2130  ptr = LALInferenceGetVariable(priorArgs, meanName);
2131  if ( ptr ) *mu = *(REAL8 *)ptr;
2132  else XLAL_ERROR_VOID(XLAL_EFAILED);
2133 
2134  ptr = LALInferenceGetVariable(priorArgs, sigmaName);
2135  if ( ptr ) *sigma = *(REAL8 *)ptr;
2136  else XLAL_ERROR_VOID(XLAL_EFAILED);
2137 
2138  ptr = LALInferenceGetVariable(priorArgs, idxName);
2139  if ( ptr ) *idx = *(UINT4 *)ptr;
2140  else XLAL_ERROR_VOID(XLAL_EFAILED);
2141 
2142  return;
2143 }
2144 
2145 /* Remove the correlated prior (for all variables) */
2146 void LALInferenceRemoveCorrelatedPrior(LALInferenceVariables *priorArgs){
2147  char corName[VARNAME_MAX];
2148  char invName[VARNAME_MAX];
2149  char meanName[VARNAME_MAX];
2150  char sigmaName[VARNAME_MAX];
2151  char idxName[VARNAME_MAX];
2152 
2153  sprintf(corName,"correlation_matrix");
2154  sprintf(invName,"inverse_correlation_matrix");
2155  sprintf(meanName, "_correlation_mean");
2156  sprintf(sigmaName, "_correlation_sigma");
2157  sprintf(idxName,"_index");
2158 
2159  LALInferenceRemoveVariable(priorArgs, corName);
2160  LALInferenceRemoveVariable(priorArgs, invName);
2161 
2162  /* remove correlated prior parameters */
2163  LALInferenceVariableItem *item = priorArgs->head;
2164  for( ; item; item = item->next ){
2165  if ( strstr(item->name, meanName) != NULL && strstr(item->name, sigmaName) != NULL && strstr(item->name, idxName) != NULL ){
2166  LALInferenceRemoveVariable(priorArgs, meanName);
2167  LALInferenceRemoveVariable(priorArgs, sigmaName);
2168  LALInferenceRemoveVariable(priorArgs, idxName);
2169  }
2170  }
2171 
2172  return;
2173 }
2174 
2175 /* Check for a correlated prior of the standard form */
2176 int LALInferenceCheckCorrelatedPrior(LALInferenceVariables *priorArgs,
2177  const char *name){
2178  char corName[VARNAME_MAX];
2179  char invName[VARNAME_MAX];
2180  char meanName[VARNAME_MAX];
2181  char sigmaName[VARNAME_MAX];
2182  char idxName[VARNAME_MAX];
2183 
2184  sprintf(corName,"correlation_matrix");
2185  sprintf(invName,"inverse_correlation_matrix");
2186  sprintf(meanName, "%s_correlation_mean", name);
2187  sprintf(sigmaName, "%s_correlation_sigma", name);
2188  sprintf(idxName,"%s_index",name);
2189 
2190  return (LALInferenceCheckVariable(priorArgs,corName) &&
2191  LALInferenceCheckVariable(priorArgs,invName) &&
2192  LALInferenceCheckVariable(priorArgs,idxName) &&
2193  LALInferenceCheckVariable(priorArgs,meanName) &&
2194  LALInferenceCheckVariable(priorArgs,sigmaName));
2195 }
2196 
2197 
2198 /* Function to add a Gaussian Mixture Model prior (without any hyperparameters)
2199  * \c name is either a single parameter name (e.g. "H0" for a 1D GMM) or a set of parameter
2200  * names separated by colons, e.g. "H0:COSIOTA" for a multivariate GMM. The means, covariance
2201  * matrices, and weights of each mode must be supplied, along with vectors containg the
2202  * minium and maximum allowed prior extent for each parameter. */
2203 void LALInferenceAddGMMPrior( LALInferenceVariables *priorArgs, const char *name,
2204  REAL8Vector ***mus, gsl_matrix ***covs,
2205  REAL8Vector **weights, REAL8Vector **minrange, REAL8Vector **maxrange ){
2206  char gmmParsName[VARNAME_MAX] = "gmm_parameter_lists"; // contains list of ':'-separated lists of GMM parameters
2207  LALStringVector *parLists = NULL;
2208 
2209  /* check if any other GMM prior already exist */
2210  if ( !LALInferenceCheckVariable( priorArgs, gmmParsName ) ){
2211  /* if not then create a new list */
2212  parLists = XLALAppendString2Vector( parLists, name );
2213  }
2214  else{
2215  /* otherwise add to existing list */
2216  parLists = *(LALStringVector **)LALInferenceGetVariable( priorArgs, gmmParsName );
2217  parLists = XLALAppendString2Vector( parLists, name );
2218  }
2219  LALInferenceAddVariable( priorArgs, gmmParsName, &parLists, LALINFERENCE_void_ptr_t, LALINFERENCE_PARAM_FIXED );
2220 
2221  char musName[VARNAME_MAX]; /* array of means for each parameter for each GMM component */
2222  char sigmasName[VARNAME_MAX]; /* array of standard deviations from each parameter for each GMM component */
2223  char weightsName[VARNAME_MAX]; /* weights for each GMM component */
2224  char corName[VARNAME_MAX]; /* correlation matrices of GMM components */
2225  char invcorName[VARNAME_MAX]; /* inverse correlation matrices of GMM components */
2226  char detName[VARNAME_MAX]; /* determinants of the GMM components */
2227  char minName[VARNAME_MAX]; /* minimum extent of prior for each parameter */
2228  char maxName[VARNAME_MAX]; /* maximum extent of prior for each parameter */
2229 
2230  sprintf(musName, "%s_gmm_sigmas", name);
2231  sprintf(sigmasName, "%s_gmm_mus", name);
2232  sprintf(weightsName, "%s_gmm_weights", name);
2233  sprintf(corName, "%s_gmm_cors", name);
2234  sprintf(invcorName, "%s_gmm_invcors", name);
2235  sprintf(detName, "%s_gmm_dets", name);
2236  sprintf(minName, "%s_gmm_min", name);
2237  sprintf(maxName, "%s_gmm_max", name);
2238 
2239  /* get number of parameters */
2240  TokenList *toks = NULL;
2241  if ( XLALCreateTokenList( &toks, name, ":" ) != XLAL_SUCCESS ){
2242  XLAL_ERROR_VOID( XLAL_EINVAL, "Could not create ':' separated token list for GMM prior" );
2243  }
2244  UINT4 npars = toks->nTokens;
2245 
2246  if ( !mus[0][0] || !covs[0][0] || !weights[0] ){
2247  XLAL_ERROR_VOID( XLAL_EINVAL, "GMM means, covariance matrices and weights must all be specified" );
2248  }
2249 
2250  /* get number of modes from weights vector */
2251  UINT4 nmodes = weights[0]->length;
2252 
2253  /* convert covariance matrix to correlation matrix (this will avoid numerical dynamic range issues when
2254  * drawing from multivariate Gaussians) */
2255  REAL8Vector **sigmas = XLALCalloc(nmodes, sizeof(REAL8Vector *)); /* get standard deviations of parameters for each mode */
2256  UINT4 j = 0, i = 0;
2257  gsl_matrix **cormat = XLALCalloc(nmodes, sizeof(gsl_matrix *)); // correlation matrices
2258  gsl_matrix **invcor = XLALCalloc(nmodes, sizeof(gsl_matrix *)); // inverse correlation matrices
2259  REAL8Vector *dets = XLALCreateREAL8Vector(nmodes); // determinants of covariance matrices
2260  for ( j = 0; j < nmodes; j++ ){
2261  sigmas[j] = XLALCreateREAL8Vector( npars );
2262  cormat[j] = gsl_matrix_calloc(npars, npars); /* initialise to all zeros */
2263  gsl_matrix *invSig = gsl_matrix_calloc(npars, npars); /* diagonal matrix containing inverse of standard deviations */
2264  gsl_matrix *tmpMat = gsl_matrix_calloc(npars, npars);
2265 
2266  if ( covs[0][j]->size1 != covs[0][j]->size2 || covs[0][j]->size1 != npars ){
2267  XLAL_ERROR_VOID( XLAL_EINVAL, "GMM covariance matrices must be square and have the correct number of parameters" );
2268  }
2269 
2270  for ( i = 0; i < npars; i++ ){
2271  sigmas[j]->data[i] = sqrt(gsl_matrix_get(covs[0][j], i, i));
2272  gsl_matrix_set(invSig, i, i, 1./sigmas[j]->data[i]);
2273  }
2274  /* get correlation coefficient matrix */
2275  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, invSig, covs[0][j], 0., tmpMat);
2276  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, tmpMat, invSig, 0., cormat[j]);
2277 
2278  gsl_matrix_free(invSig);
2279 
2280  /* get inverse of correlation matrix */
2281  XLAL_CALLGSL( gsl_matrix_memcpy( tmpMat, cormat[j] ) ); // copy of correlation matrix
2282 
2283  invcor[j] = gsl_matrix_alloc( npars, npars );
2284  gsl_permutation *p = gsl_permutation_alloc( npars );
2285  INT4 s;
2286 
2287  /* check correlation matrix is positive definite */
2288  if( !LALInferenceCheckPositiveDefinite( tmpMat, tmpMat->size1 ) ){
2289  XLAL_ERROR_VOID( XLAL_EFUNC | XLAL_EINVAL, "Error... matrix is not positive definite!" );
2290  }
2291 
2292  /* invert correlation matrix */
2293  XLAL_CALLGSL( gsl_linalg_LU_decomp( tmpMat, p, &s ) );
2294  XLAL_CALLGSL( gsl_linalg_LU_invert( tmpMat, p, invcor[j] ) );
2295 
2296  /* get determinant of covariance matrix */
2297  dets->data[j] = gsl_linalg_LU_det( tmpMat, s );
2298 
2299  XLAL_CALLGSL( gsl_permutation_free( p ) );
2300  XLAL_CALLGSL( gsl_matrix_free( tmpMat ) );
2301  }
2302 
2303  /* add values */
2304  LALInferenceAddVariable( priorArgs, musName, mus, LALINFERENCE_void_ptr_t, LALINFERENCE_PARAM_FIXED );
2305  LALInferenceAddVariable( priorArgs, sigmasName, &sigmas, LALINFERENCE_void_ptr_t, LALINFERENCE_PARAM_FIXED );
2306  LALInferenceAddVariable( priorArgs, corName, &cormat, LALINFERENCE_void_ptr_t, LALINFERENCE_PARAM_FIXED );
2307  LALInferenceAddVariable( priorArgs, invcorName, &invcor, LALINFERENCE_void_ptr_t, LALINFERENCE_PARAM_FIXED );
2308  LALInferenceAddVariable( priorArgs, detName, &dets, LALINFERENCE_REAL8Vector_t, LALINFERENCE_PARAM_FIXED );
2309 
2310  /* make sure weights are normalised to 1 */
2311  REAL8 weightsum = 0.;
2312  for ( i = 0; i < weights[0]->length; i++ ){ weightsum += weights[0]->data[i]; }
2313  for ( i = 0; i < weights[0]->length; i++ ){ weights[0]->data[i] /= weightsum; }
2314 
2315  LALInferenceAddVariable( priorArgs, weightsName, weights, LALINFERENCE_REAL8Vector_t, LALINFERENCE_PARAM_FIXED );
2316 
2317  /* check if ranges are set and if not set to +/- infinity */
2318  if ( !minrange[0] ){
2319  minrange[0] = XLALCreateREAL8Vector( npars );
2320  for ( i = 0; i < npars; i++ ){ minrange[0]->data[i] = -INFINITY; } /* set lower bound to -infinity */
2321  }
2322  if ( !maxrange[0] ){
2323  maxrange[0] = XLALCreateREAL8Vector( npars );
2324  for ( i = 0; i < npars; i++ ){ maxrange[0]->data[i] = INFINITY; } /* set upper bound to infinity */
2325  }
2326 
2327  for ( i = 0; i < npars; i++ ){
2328  if ( minrange[0]->data[i] > maxrange[0]->data[i] ){
2329  XLAL_ERROR_VOID( XLAL_EINVAL, "Bounds for GMM are wrong" );
2330  }
2331  }
2332 
2333  LALInferenceAddVariable( priorArgs, minName, minrange, LALINFERENCE_REAL8Vector_t, LALINFERENCE_PARAM_FIXED );
2334  LALInferenceAddVariable( priorArgs, maxName, maxrange, LALINFERENCE_REAL8Vector_t, LALINFERENCE_PARAM_FIXED );
2335  return;
2336 }
2337 
2338 /* get the standard deviations, correlation matrices, inverse correlation matrices, means, weights and bounds
2339  * of the GMM prior, the index of the parameter given in name, and the full (i.e. including all parameters)
2340  * name of the GMM prior */
2341 void LALInferenceGetGMMPrior( LALInferenceVariables *priorArgs, const char *name,
2342  REAL8Vector ***mus, REAL8Vector ***sigmas, gsl_matrix ***cors, gsl_matrix ***invcors,
2343  REAL8Vector **weights, REAL8Vector **minrange, REAL8Vector **maxrange,
2344  REAL8Vector **dets, UINT4 *idx, CHAR **fullname ){
2345  /* find list of GMM parameters that this parameter lives in */
2346  char gmmParsName[VARNAME_MAX] = "gmm_parameter_lists"; // contains list of ':'-separated lists of GMM parameters
2347 
2348  /* check if the given variable is in any of the lists */
2349  LALStringVector *parLists = *(LALStringVector **)LALInferenceGetVariable( priorArgs, gmmParsName );
2350  UINT4 numlists = parLists->length;
2351 
2352  UINT4 found = 0, i = 0, j = 0;
2353  for ( i = 0; i < numlists; i++ ){
2354  TokenList *toks = NULL;
2355  if ( XLALCreateTokenList( &toks, parLists->data[i], ":" ) != XLAL_SUCCESS ){ XLAL_ERROR_VOID(XLAL_EFAILED); }
2356  for ( j = 0; j < toks->nTokens; j++ ){
2357  if ( !strcmp(toks->tokens[j], name) ){
2358  found = 1;
2359  *idx = j; // set index of parameter
2360  break;
2361  }
2362  }
2363  XLALDestroyTokenList( toks );
2364  if ( found ){ break; }
2365  }
2366  if ( !found ){ XLAL_ERROR_VOID(XLAL_EFAILED); } // parameter could not be found
2367 
2368  *fullname = parLists->data[i];
2369 
2370  char musName[VARNAME_MAX];
2371  char sigmasName[VARNAME_MAX];
2372  char weightsName[VARNAME_MAX];
2373  char corName[VARNAME_MAX];
2374  char invcorName[VARNAME_MAX];
2375  char detName[VARNAME_MAX];
2376  char minName[VARNAME_MAX];
2377  char maxName[VARNAME_MAX];
2378  void *ptr = NULL;
2379 
2380  sprintf(musName, "%s_gmm_sigmas", parLists->data[i]);
2381  sprintf(sigmasName, "%s_gmm_mus", parLists->data[i]);
2382  sprintf(weightsName, "%s_gmm_weights", parLists->data[i]);
2383  sprintf(corName, "%s_gmm_cors", parLists->data[i]);
2384  sprintf(invcorName, "%s_gmm_invcors", parLists->data[i]);
2385  sprintf(detName, "%s_gmm_dets", parLists->data[i]);
2386  sprintf(minName, "%s_gmm_min", parLists->data[i]);
2387  sprintf(maxName, "%s_gmm_max",parLists->data[i]);
2388 
2389  ptr = LALInferenceGetVariable(priorArgs, musName);
2390  if ( ptr ){ *mus = *(REAL8Vector ***)ptr; }
2391  else{ XLAL_ERROR_VOID(XLAL_EFAILED); }
2392 
2393  ptr = LALInferenceGetVariable(priorArgs, sigmasName);
2394  if ( ptr ){ *sigmas = *(REAL8Vector ***)ptr; }
2395  else{ XLAL_ERROR_VOID(XLAL_EFAILED); }
2396 
2397  ptr = LALInferenceGetVariable(priorArgs, weightsName);
2398  if ( ptr ){ *weights = *(REAL8Vector **)ptr; }
2399  else{ XLAL_ERROR_VOID(XLAL_EFAILED); }
2400 
2401  ptr = LALInferenceGetVariable(priorArgs, corName);
2402  if ( ptr ){ *cors = *(gsl_matrix ***)ptr; }
2403  else{ XLAL_ERROR_VOID(XLAL_EFAILED); }
2404 
2405  ptr = LALInferenceGetVariable(priorArgs, invcorName);
2406  if ( ptr ){ *invcors = *(gsl_matrix ***)ptr; }
2407  else{ XLAL_ERROR_VOID(XLAL_EFAILED); }
2408 
2409  ptr = LALInferenceGetVariable(priorArgs, detName);
2410  if ( ptr ){ *dets = *(REAL8Vector **)ptr; }
2411  else{ XLAL_ERROR_VOID(XLAL_EFAILED); }
2412 
2413  ptr = LALInferenceGetVariable(priorArgs, minName);
2414  if ( ptr ){ *minrange = *(REAL8Vector **)ptr; }
2415  else{ XLAL_ERROR_VOID(XLAL_EFAILED); }
2416 
2417  ptr = LALInferenceGetVariable(priorArgs, maxName);
2418  if ( ptr ){ *maxrange = *(REAL8Vector **)ptr; }
2419  else{ XLAL_ERROR_VOID(XLAL_EFAILED); }
2420  return;
2421 }
2422 
2423 /* check for GMM prior */
2424 int LALInferenceCheckGMMPrior(LALInferenceVariables *priorArgs, const char *name){
2425  char gmmParsName[VARNAME_MAX] = "gmm_parameter_lists"; // contains list of ':'-separated lists of GMM parameters
2426 
2427  /* check if there is a list of GMM parameters */
2428  if ( !LALInferenceCheckVariable(priorArgs, gmmParsName) ){ return 0; }
2429 
2430  /* check if the given variable is in any of the lists */
2431  LALStringVector *parLists = *(LALStringVector **)LALInferenceGetVariable( priorArgs, gmmParsName );
2432  UINT4 numlists = parLists->length;
2433 
2434  if ( numlists == 0 ){ return 0; }
2435 
2436  UINT4 found = 0, i = 0, j = 0;
2437  for ( i = 0; i < numlists; i++ ){
2438  TokenList *toks = NULL;
2439  if ( XLALCreateTokenList( &toks, parLists->data[i], ":" ) != XLAL_SUCCESS ){ return 0; }
2440  for ( j = 0; j < toks->nTokens; j++ ){
2441  if ( !strcmp(toks->tokens[j], name) ){
2442  found = 1;
2443  break;
2444  }
2445  }
2446  XLALDestroyTokenList( toks );
2447  if ( found ){ break; }
2448  }
2449 
2450  if ( !found ){ return 0; } // parameter could not be found
2451 
2452  /* check for required items for the GMM containing the found parameter */
2453  char musName[VARNAME_MAX];
2454  char sigmasName[VARNAME_MAX];
2455  char weightsName[VARNAME_MAX];
2456  char corName[VARNAME_MAX];
2457  char invcorName[VARNAME_MAX];
2458  char detName[VARNAME_MAX];
2459  char minName[VARNAME_MAX];
2460  char maxName[VARNAME_MAX];
2461 
2462  sprintf(musName, "%s_gmm_sigmas", parLists->data[i]);
2463  sprintf(sigmasName, "%s_gmm_mus", parLists->data[i]);
2464  sprintf(weightsName, "%s_gmm_weights", parLists->data[i]);
2465  sprintf(corName, "%s_gmm_cors", parLists->data[i]);
2466  sprintf(invcorName, "%s_gmm_invcors", parLists->data[i]);
2467  sprintf(detName, "%s_gmm_dets", parLists->data[i]);
2468  sprintf(minName, "%s_gmm_min", parLists->data[i]);
2469  sprintf(maxName, "%s_gmm_max", parLists->data[i]);
2470 
2471  return (LALInferenceCheckVariable(priorArgs, musName) &&
2472  LALInferenceCheckVariable(priorArgs, sigmasName) &&
2473  LALInferenceCheckVariable(priorArgs, weightsName) &&
2474  LALInferenceCheckVariable(priorArgs, corName) &&
2475  LALInferenceCheckVariable(priorArgs, invcorName) &&
2476  LALInferenceCheckVariable(priorArgs, minName) &&
2477  LALInferenceCheckVariable(priorArgs, maxName) &&
2478  LALInferenceCheckVariable(priorArgs, detName));
2479 }
2480 
2481 
2482 /* remove GMM Prior (either using the full name [i.e. comma separated list of GMM parameters], or
2483  * a single parameter names, from the GMM prior. Note that the corresponding entry in the "gmm_parameter_lists"
2484  * will be set to a string termination character, or if it just contains one entry it will be removed. */
2485 void LALInferenceRemoveGMMPrior( LALInferenceVariables *priorArgs, const char *name ){
2486  char gmmParsName[VARNAME_MAX] = "gmm_parameter_lists"; // contains list of ':'-separated lists of GMM parameters
2487 
2488  /* check if the given variable is in any of the lists */
2489  LALStringVector *parLists = *(LALStringVector **)LALInferenceGetVariable( priorArgs, gmmParsName );
2490  UINT4 numlists = parLists->length;
2491 
2492  UINT4 found = 0, i = 0, j = 0;
2493  for ( i = 0; i < numlists; i++ ){
2494  if ( !strcmp(name, parLists->data[i]) ){
2495  found = 1;
2496  break;
2497  }
2498  TokenList *toks = NULL;
2499  if ( XLALCreateTokenList( &toks, parLists->data[i], ":" ) != XLAL_SUCCESS ){ XLAL_ERROR_VOID(XLAL_EFAILED); }
2500  for ( j = 0; j < toks->nTokens; j++ ){
2501  if ( !strcmp(toks->tokens[j], name) ){
2502  found = 1;
2503  break;
2504  }
2505  }
2506  XLALDestroyTokenList( toks );
2507  if ( found ){ break; }
2508  }
2509 
2510  if ( !found ){ return; } // parameter doesn't seem to be present anyway!
2511 
2512  char musName[VARNAME_MAX];
2513  char sigmasName[VARNAME_MAX];
2514  char weightsName[VARNAME_MAX];
2515  char corName[VARNAME_MAX];
2516  char invcorName[VARNAME_MAX];
2517  char detName[VARNAME_MAX];
2518  char minName[VARNAME_MAX];
2519  char maxName[VARNAME_MAX];
2520 
2521  sprintf(musName, "%s_gmm_sigmas", parLists->data[i]);
2522  sprintf(sigmasName, "%s_gmm_mus", parLists->data[i]);
2523  sprintf(weightsName, "%s_gmm_weights", parLists->data[i]);
2524  sprintf(corName, "%s_gmm_cors", parLists->data[i]);
2525  sprintf(invcorName, "%s_gmm_invcors", parLists->data[i]);
2526  sprintf(detName, "%s_gmm_dets", parLists->data[i]);
2527  sprintf(minName, "%s_gmm_min", parLists->data[i]);
2528  sprintf(maxName, "%s_gmm_max", parLists->data[i]);
2529 
2530  LALInferenceRemoveVariable(priorArgs, musName);
2531  LALInferenceRemoveVariable(priorArgs, sigmasName);
2532  LALInferenceRemoveVariable(priorArgs, weightsName);
2533  LALInferenceRemoveVariable(priorArgs, corName);
2534  LALInferenceRemoveVariable(priorArgs, invcorName);
2535  LALInferenceRemoveVariable(priorArgs, detName);
2536  LALInferenceRemoveVariable(priorArgs, minName);
2537  LALInferenceRemoveVariable(priorArgs, maxName);
2538 
2539  if ( numlists == 1 ){ /* just remove list */
2540  LALInferenceRemoveVariable(priorArgs, gmmParsName);
2541  }
2542  else{ /* set found entry in "gmm_parameter_lists" to string termination character */
2543  parLists->data[i][0] = '\0';
2544  }
2545 
2546  return;
2547 }
2548 
2549 
2550 /* Check for a Fermi-Dirac Prior */
2551 int LALInferenceCheckFermiDiracPrior(LALInferenceVariables *priorArgs, const char *name)
2552 {
2553  char rName[VARNAME_MAX];
2554  char sigmaName[VARNAME_MAX];
2555  sprintf(rName,"%s_fermi_r",name);
2556  sprintf(sigmaName,"%s_fermi_sigma",name);
2557  return (LALInferenceCheckVariable(priorArgs,rName) && LALInferenceCheckVariable(priorArgs,sigmaName));
2558 }
2559 
2560 /* Function to add the r and sigma values for the prior onto the priorArgs */
2561 void LALInferenceAddFermiDiracPrior( LALInferenceVariables *priorArgs,
2562  const char *name, REAL8 *sigma, REAL8 *r,
2563  LALInferenceVariableType type ){
2564  char rName[VARNAME_MAX];
2565  char sigmaName[VARNAME_MAX];
2566 
2567  sprintf(rName,"%s_fermi_r",name);
2568  sprintf(sigmaName,"%s_fermi_sigma",name);
2569 
2570  LALInferenceAddVariable(priorArgs,rName,r,type,LALINFERENCE_PARAM_FIXED);
2571  LALInferenceAddVariable(priorArgs,sigmaName,sigma,type,LALINFERENCE_PARAM_FIXED);
2572  return;
2573 }
2574 
2575 /* Function to remove the r and sigma values for the prior onto the priorArgs */
2576 void LALInferenceRemoveFermiDiracPrior(LALInferenceVariables *priorArgs, const char *name){
2577  char rName[VARNAME_MAX];
2578  char sigmaName[VARNAME_MAX];
2579 
2580  sprintf(rName,"%s_fermi_r",name);
2581  sprintf(sigmaName,"%s_fermi_sigma",name);
2582 
2583  LALInferenceRemoveVariable(priorArgs, rName);
2584  LALInferenceRemoveVariable(priorArgs, sigmaName);
2585  return;
2586 }
2587 
2588 /* Get the r and sigma values of the prior from the priorArgs list, given a name */
2589 void LALInferenceGetFermiDiracPrior(LALInferenceVariables *priorArgs,
2590  const char *name, REAL8 *sigma, REAL8 *r)
2591 {
2592  char rName[VARNAME_MAX];
2593  char sigmaName[VARNAME_MAX];
2594  void *ptr=NULL;
2595 
2596  sprintf(rName,"%s_fermi_r",name);
2597  sprintf(sigmaName,"%s_fermi_sigma",name);
2598 
2599  ptr = LALInferenceGetVariable(priorArgs, rName);
2600  if ( ptr ) *r = *(REAL8*)ptr;
2601  else XLAL_ERROR_VOID(XLAL_EFAILED);
2602 
2603  ptr = LALInferenceGetVariable(priorArgs, sigmaName);
2604  if ( ptr ) *sigma = *(REAL8*)ptr;
2605  else XLAL_ERROR_VOID(XLAL_EFAILED);
2606 
2607  return;
2608 }
2609 
2610 
2611 /* Check for a prior uniform in the log */
2612 int LALInferenceCheckLogUniformPrior(LALInferenceVariables *priorArgs,
2613  const char *name)
2614 {
2615  char xminName[VARNAME_MAX];
2616  char xmaxName[VARNAME_MAX];
2617 
2618  sprintf(xminName, "%s_loguniform_xmin", name);
2619  sprintf(xmaxName, "%s_loguniform_xmax", name);
2620 
2621  return (LALInferenceCheckVariable(priorArgs, xminName) &&
2622  LALInferenceCheckVariable(priorArgs, xmaxName));
2623 }
2624 
2625 /* Add domain boundary values for the prior onto the priorArgs */
2626 void LALInferenceAddLogUniformPrior(LALInferenceVariables *priorArgs,
2627  const char *name, REAL8 *xmin, REAL8 *xmax,
2628  LALInferenceVariableType type )
2629 {
2630  if (*xmin >= *xmax || *xmin < 0. ){
2631  XLAL_ERROR_VOID(XLAL_EINVAL, "Minimum must be less than maximum (and minumum must be greater than zero), but %f >= %f.", *xmin, *xmax);
2632  }
2633 
2634  char xminName[VARNAME_MAX];
2635  char xmaxName[VARNAME_MAX];
2636 
2637  sprintf(xminName, "%s_loguniform_xmin", name);
2638  sprintf(xmaxName, "%s_loguniform_xmax", name);
2639 
2640  LALInferenceAddVariable(priorArgs, xminName, xmin, type,
2641  LALINFERENCE_PARAM_FIXED);
2642  LALInferenceAddVariable(priorArgs, xmaxName, xmax, type,
2643  LALINFERENCE_PARAM_FIXED);
2644  return;
2645 }
2646 
2647 /* Remove the domain boundary values for the prior onto the priorArgs */
2648 void LALInferenceRemoveLogUniformPrior(LALInferenceVariables *priorArgs,
2649  const char *name)
2650 {
2651  char xminName[VARNAME_MAX];
2652  char xmaxName[VARNAME_MAX];
2653 
2654  sprintf(xminName, "%s_loguniform_xmin", name);
2655  sprintf(xmaxName, "%s_loguniform_xmax", name);
2656 
2657  LALInferenceRemoveVariable(priorArgs, xminName);
2658  LALInferenceRemoveVariable(priorArgs, xmaxName);
2659  return;
2660 }
2661 
2662 /* Get domain boundary values of the prior from priorArgs list, given a name */
2663 void LALInferenceGetLogUniformPrior(LALInferenceVariables *priorArgs,
2664  const char *name, REAL8 *xmin, REAL8 *xmax)
2665 {
2666  char xminName[VARNAME_MAX];
2667  char xmaxName[VARNAME_MAX];
2668  void *ptr=NULL;
2669 
2670  sprintf(xminName, "%s_loguniform_xmin", name);
2671  sprintf(xmaxName, "%s_loguniform_xmax", name);
2672 
2673  ptr = LALInferenceGetVariable(priorArgs, xminName);
2674  if ( ptr ) {
2675  *xmin = *(REAL8*)ptr;
2676  } else {
2677  XLAL_ERROR_VOID(XLAL_EFAILED);
2678  }
2679 
2680  ptr = LALInferenceGetVariable(priorArgs, xmaxName);
2681  if ( ptr ) {
2682  *xmax = *(REAL8*)ptr;
2683  } else {
2684  XLAL_ERROR_VOID(XLAL_EFAILED);
2685  }
2686 
2687  return;
2688 }
2689 
2690 void LALInferenceDrawFromPrior( LALInferenceVariables *output,
2691  LALInferenceVariables *priorArgs,
2692  gsl_rng *rdm) {
2693  if (output == NULL || priorArgs == NULL || rdm == NULL)
2694  XLAL_ERROR_VOID(XLAL_EFAULT, "Null arguments received.");
2695 
2696  LALInferenceVariableItem *item = output->head;
2697 
2698  /* check if using a k-D tree as the prior */
2699  if( LALInferenceCheckVariable( priorArgs, "kDTreePrior" ) ){
2700  LALInferenceKDTree *tree =
2701  *(LALInferenceKDTree **)LALInferenceGetVariable(priorArgs, "kDTreePrior");
2702 
2703  /* get parameter template */
2704  LALInferenceVariables *templt =
2705  *(LALInferenceVariables **)LALInferenceGetVariable(priorArgs,
2706  "kDTreePriorTemplate");
2707 
2708  UINT4 Ncell = 8; /* number of points in a prior cell - i.e. controls
2709  how fine or coarse the prior looks (default to 8) */
2710 
2711  if( LALInferenceCheckVariable( priorArgs, "kDTreePriorNcell" ) )
2712  Ncell = *(UINT4 *)LALInferenceGetVariable( priorArgs,"kDTreePriorNcell");
2713 
2714  /* draw all points from the prior distribution */
2715  REAL8 *proposedPt = XLALCalloc(tree->dim, sizeof(REAL8));
2716 
2717  /* A randomly-chosen point from those in the tree. */
2718  LALInferenceKDDrawEigenFrame(rdm, tree, proposedPt, Ncell);
2719  LALInferenceKDREAL8ToVariables(output, proposedPt, templt);
2720  }
2721  else{
2722  for(;item;item=item->next){
2723  if(item->vary==LALINFERENCE_PARAM_CIRCULAR ||
2724  item->vary==LALINFERENCE_PARAM_LINEAR)
2725  LALInferenceDrawNameFromPrior( output, priorArgs, item->name,
2726  item->type, rdm );
2727  }
2728 
2729  /* remove multivariate deviates value if set */
2730  if ( LALInferenceCheckVariable( priorArgs, "multivariate_deviates" ) )
2731  LALInferenceRemoveVariable( priorArgs, "multivariate_deviates" );
2732 
2733  }
2734 }
2735 
2736 void LALInferenceDrawNameFromPrior( LALInferenceVariables *output,
2737  LALInferenceVariables *priorArgs,
2738  char *name, LALInferenceVariableType type,
2739  gsl_rng *rdm) {
2740  if (output == NULL || priorArgs == NULL || name == NULL || rdm == NULL)
2741  XLAL_ERROR_VOID(XLAL_EFAULT, "Null arguments received.");
2742 
2743  REAL8 tmp = 0.;
2744 
2745  /* test for a Gaussian prior */
2746  if( LALInferenceCheckGaussianPrior( priorArgs, name ) ){
2747  REAL8 mu = 0., sigma = 0.;
2748 
2749  LALInferenceGetGaussianPrior( priorArgs, name, (void *)&mu, (void *)&sigma );
2750  tmp = mu + gsl_ran_gaussian(rdm, (double)sigma);
2751  }
2752  /* test for uniform prior */
2753  else if( LALInferenceCheckMinMaxPrior( priorArgs, name ) ){
2754  REAL8 min = 0., max = 0.;
2755 
2756  LALInferenceGetMinMaxPrior(priorArgs, name, &min, &max);
2757  tmp = min + (max-min)*gsl_rng_uniform( rdm );
2758  }
2759  /* test for a Fermi-Dirac prior */
2760  else if( LALInferenceCheckFermiDiracPrior( priorArgs, name ) ){
2761  REAL8 r = 0., sigma = 0., cp;
2762 
2763  LALInferenceGetFermiDiracPrior(priorArgs, name, &sigma, &r);
2764 
2765  /* use the inverse sampling transform to draw a new sample */
2766  do { /* numerical issues mean that the analytic solution to this equation can go negative, so make sure that is not the case */
2767  cp = gsl_rng_uniform( rdm ); /* draw a point uniformly between 0 and 1 */
2768  tmp = log(-exp(-r) + pow(1. + exp(r), -cp) + exp(-r)*pow(1. + exp(r), -cp));
2769  tmp *= -sigma;
2770  } while ( tmp < 0. );
2771  }
2772  /* test for a prior uniform in the log */
2773  else if( LALInferenceCheckLogUniformPrior( priorArgs, name ) ){
2774  REAL8 xmin = 0., xmax = 0., cp;
2775 
2776  LALInferenceGetLogUniformPrior(priorArgs, name, &xmin, &xmax);
2777 
2778  if( xmin <= 0 ) {
2779  XLAL_ERROR_VOID(XLAL_EDOM, "Log-uniform min value not positive.");
2780  } else if( xmax < xmin ) {
2781  XLAL_ERROR_VOID(XLAL_EDOM, "Log-uniform min value greater than max.");
2782  }
2783 
2784  /* use the inverse sampling transform to draw a new sample */
2785  cp = gsl_rng_uniform( rdm ); /* draw a point uniformly between 0 and 1 */
2786 
2787  /* the percentage-point function (PPF, or inverse CDF) for PDF~1/x is:
2788  \f$x = (\frac{x_{\rm max}}{x_{\rm min}})^{\rm cdf} x_{\rm min}\f$ */
2789  tmp = xmin * pow(xmax/xmin, cp);
2790  }
2791  /* test for a prior drawn from correlated values */
2792  else if( LALInferenceCheckCorrelatedPrior( priorArgs, name ) ){
2793  gsl_matrix *cor = NULL, *invcor = NULL;
2794  REAL8 mu = 0, sigma = 0.;
2795  UINT4 idx = 0, dims = 0;
2796  REAL4Vector *tmps = NULL;
2797 
2798  LALInferenceGetCorrelatedPrior( priorArgs, name, &cor, &invcor, &mu, &sigma, &idx );
2799  dims = cor->size1;
2800 
2801  /* to avoid unnecessary repetition the multivariate deviates are be
2802  added as a new variable that can be extracted during multiple calls.
2803  This will be later removed by the LALInferenceDrawFromPrior function. */
2804  if ( LALInferenceCheckVariable( priorArgs, "multivariate_deviates" ) ){
2805  tmps = *(REAL4Vector **)LALInferenceGetVariable(priorArgs,
2806  "multivariate_deviates");
2807  }
2808  else{
2809  RandomParams *randParam;
2810  UINT4 randomseed = gsl_rng_get(rdm);
2811 
2812  /* check matrix for positive definiteness */
2813  if( !LALInferenceCheckPositiveDefinite( cor, dims ) ){
2814  XLAL_ERROR_VOID(XLAL_EFUNC | XLAL_EINVAL, "Matrix is not positive-definite!");
2815  }
2816 
2817  /* draw values from the multivariate Gaussian described by the correlation matrix */
2818  tmps = XLALCreateREAL4Vector( dims );
2819  randParam = XLALCreateRandomParams( randomseed );
2820  XLALMultiNormalDeviates( tmps, cor, dims, randParam );
2821 
2822  LALInferenceAddVariable( priorArgs, "multivariate_deviates", &tmps,
2823  LALINFERENCE_REAL8Vector_t,
2824  LALINFERENCE_PARAM_FIXED );
2825  XLALDestroyRandomParams( randParam );
2826  }
2827 
2828  /* set random number for given parameter index (converted to a draw from the covariance matrix) */
2829  tmp = mu + sigma*tmps->data[idx];
2830 
2831  /* free tmps */
2832  if ( !LALInferenceCheckVariable( priorArgs, "multivariate_deviates" ) )
2833  XLALDestroyREAL4Vector( tmps );
2834  }
2835  /* test if a prior drawn from a Gaussian Mixture Model */
2836  else if ( LALInferenceCheckGMMPrior( priorArgs, name ) ){
2837  gsl_matrix **cor, UNUSED **invcor;
2838  REAL8Vector **mus, **sigmas, *weights = NULL, *minrange = NULL, *maxrange = NULL, UNUSED *dets = NULL;
2839  REAL8 cp = 0., cumweights = 0.;
2840  REAL4Vector *tmps = NULL;
2841  UINT4 idx = 0, dims = 0;
2842  CHAR *fullname = NULL;
2843 
2844  LALInferenceGetGMMPrior( priorArgs, name, &mus, &sigmas, &cor, &invcor, &weights, &minrange, &maxrange, &dets, &idx, &fullname );
2845 
2846  dims = cor[0]->size1;
2847 
2848  /* set mode to use and multivariate deviates for it */
2849  CHAR gmmmvd[VARNAME_MAX], gmmmode[VARNAME_MAX];
2850  UINT4 thismode = 0;
2851  sprintf(gmmmvd, "%s_gmm_mvd", fullname);
2852  sprintf(gmmmode, "%s_gmm_mode", fullname);
2853  if ( LALInferenceCheckVariable( priorArgs, gmmmvd ) ){ /* get values if already set */
2854  tmps = *(REAL4Vector **)LALInferenceGetVariable(priorArgs, gmmmvd);
2855  thismode = LALInferenceGetUINT4Variable( priorArgs, gmmmode);
2856  }
2857  else{ /* otherwise create array */
2858  tmps = XLALCreateREAL4Vector( dims );
2859 
2860  /* get GMM mode to use */
2861  cp = gsl_rng_uniform( rdm );
2862  cumweights = weights->data[0];
2863  /* get index of mode */
2864  while( cp > cumweights ){
2865  thismode++;
2866  cumweights += weights->data[thismode];
2867  }
2868 
2869  RandomParams *randParam;
2870  UINT4 randomseed = gsl_rng_get(rdm);
2871  randParam = XLALCreateRandomParams( randomseed );
2872  XLALMultiNormalDeviates( tmps, cor[thismode], dims, randParam );
2873  XLALDestroyRandomParams( randParam );
2874 
2875  /* add to priorArgs */
2876  LALInferenceAddVariable( priorArgs, gmmmvd, &tmps, LALINFERENCE_void_ptr_t, LALINFERENCE_PARAM_FIXED );
2877 
2878  /* add the mode to the prior args */
2879  LALInferenceAddVariable( priorArgs, gmmmode, &thismode, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED );
2880  }
2881 
2882  /* set random number for given parameter index (converted to a draw from the covariance matrix) */
2883  tmp = mus[thismode]->data[idx] + sigmas[thismode]->data[idx]*tmps->data[idx];
2884 
2885  /* get number of parameter from GMM already passed */
2886  CHAR nparsdonename[VARNAME_MAX];
2887  UINT4 npars = 1;
2888  sprintf(nparsdonename, "%s_gmm_npars", fullname);
2889  if ( LALInferenceCheckVariable( priorArgs, nparsdonename ) ){
2890  UINT4 nparsdone = LALInferenceGetUINT4Variable( priorArgs, nparsdonename );
2891  npars = nparsdone;
2892  npars++;
2893  LALInferenceRemoveVariable( priorArgs, nparsdonename );
2894  }
2895 
2896  LALInferenceAddVariable( priorArgs, nparsdonename, &npars, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED );
2897 
2898  /* remove GMM multivariate deviates if all parameter in the GMM have been set */
2899  if ( npars == dims ){
2900  LALInferenceRemoveVariable( priorArgs, gmmmvd );
2901  LALInferenceRemoveVariable( priorArgs, gmmmode );
2902  LALInferenceRemoveVariable( priorArgs, nparsdonename );
2903  }
2904  }
2905  /* not a recognised prior type */
2906  else{
2907  return;
2908  }
2909 
2910  switch ( type ){
2911  case LALINFERENCE_REAL4_t:
2912  {
2913  REAL4 val = (REAL4)tmp;
2914  LALInferenceSetVariable(output, name, &val);
2915  break;
2916  }
2917  case LALINFERENCE_REAL8_t:
2918  {
2919  REAL8 val = tmp;
2920  LALInferenceSetVariable(output, name, &val);
2921  break;
2922  }
2923  case LALINFERENCE_INT4_t:
2924  {
2925  INT4 val = (INT4)tmp;
2926  LALInferenceSetVariable(output, name, &val);
2927  break;
2928  }
2929  case LALINFERENCE_INT8_t:
2930  {
2931  INT8 val = (INT8)tmp;
2932  LALInferenceSetVariable(output, name, &val);
2933  break;
2934  }
2935  //LALInferenceDrawFromPrior() does not handle gsl matrices properly at the moment
2936  /*case LALINFERENCE_gslMatrix_t:
2937  {
2938  REAL8 val = tmp;
2939  LALInferenceSetVariable(output, name, &val);
2940  break;
2941  }*/
2942  default:
2943  XLALPrintWarning("Trying to randomise a non-numeric parameter %s!\n",name);
2944  break;
2945  }
2946 }
2947 
2948 
2949 /* Switch reads true if parameters lie within Malmquist prior */
2950 UINT4 within_malmquist(LALInferenceRunState *runState, LALInferenceVariables *params, LALInferenceModel *model) {
2951  UINT4 i=0, nifo=0;
2952 
2953  LALInferenceIFOData *ifo = runState->data;
2954  while (ifo != NULL) {
2955  nifo++;
2956  ifo = ifo->next;
2957  }
2958 
2959  LALInferenceNetworkSNR(params, runState->data, model);
2960  REAL8 loudest_snr=0.0, second_loudest_snr=0.0;
2961  for (i=0; i<nifo; i++) {
2962  if (model->ifo_SNRs[i] > second_loudest_snr) {
2963  if (model->ifo_SNRs[i] > loudest_snr) {
2964  second_loudest_snr = loudest_snr;
2965  loudest_snr = model->ifo_SNRs[i];
2966  } else {
2967  second_loudest_snr = model->ifo_SNRs[i];
2968  }
2969  }
2970  }
2971 
2972  REAL8 malmquist_loudest = (*(REAL8 *)LALInferenceGetVariable(runState->priorArgs,"malmquist_loudest_snr"));
2973  REAL8 malmquist_second_loudest = (*(REAL8 *)LALInferenceGetVariable(runState->priorArgs,"malmquist_second_loudest_snr"));
2974  REAL8 malmquist_network = (*(REAL8 *)LALInferenceGetVariable(runState->priorArgs,"malmquist_network_snr"));
2975 
2976  if (loudest_snr < malmquist_loudest
2977  || second_loudest_snr < malmquist_second_loudest
2978  || model->SNR < malmquist_network)
2979  return(0);
2980  else
2981  return(1);
2982 }
2983 
2984 
2985 REAL8 LALInferenceAnalyticNullPrior(LALInferenceRunState UNUSED *runState, LALInferenceVariables *params, LALInferenceModel UNUSED *model) {
2986  REAL8 logPrior=0.0;
2987  REAL8 logmc=0.0,mc=0.0;
2988  REAL8 m1=0.0,m2=0.0,q=0.0,eta=0.0;
2989 
2990  if(LALInferenceCheckVariable(params,"eta")||LALInferenceCheckVariable(params,"q")) {
2991  if(LALInferenceCheckVariable(params,"logmc")) {
2992  logmc=*(REAL8 *)LALInferenceGetVariable(params,"logmc");
2993  if(LALInferenceCheckVariable(params,"q")) {
2994  q=*(REAL8 *)LALInferenceGetVariable(params,"q");
2995  LALInferenceMcQ2Masses(exp(logmc),q,&m1,&m2);
2996  logPrior+=log(m1*m1);
2997  } else {
2998  eta=*(REAL8 *)LALInferenceGetVariable(params,"eta");
2999  LALInferenceMcEta2Masses(exp(logmc),eta,&m1,&m2);
3000  logPrior+=log(((m1+m2)*(m1+m2)*(m1+m2))/(m1-m2));
3001  }
3002  /*careful using LALInferenceMcEta2Masses, it returns m1>=m2*/
3003  } else if(LALInferenceCheckVariable(params,"chirpmass")) {
3004  mc=*(REAL8 *)LALInferenceGetVariable(params,"chirpmass");
3005  if(LALInferenceCheckVariable(params,"q")) {
3006  q=*(REAL8 *)LALInferenceGetVariable(params,"q");
3007  LALInferenceMcQ2Masses(mc,q,&m1,&m2);
3008  logPrior+=log(m1*m1/mc);
3009  } else {
3010  eta=*(REAL8 *)LALInferenceGetVariable(params,"eta");
3011  LALInferenceMcEta2Masses(mc,eta,&m1,&m2);
3012  logPrior+=log(((m1+m2)*(m1+m2))/((m1-m2)*pow(eta,3.0/5.0)));
3013  }
3014  }
3015  }
3016  logPrior+=LALInferenceFlatBoundedPrior(runState, params);
3017  return(logPrior);
3018 }
3019 
3020 UINT4 LALInferenceAnalyticCubeToPrior(LALInferenceRunState *runState, LALInferenceVariables *params, UNUSED LALInferenceModel *model, double *Cube, void *context) {
3021  int i = 0;
3022  REAL8 min=-INFINITY,max=INFINITY;
3023  REAL8 m1=1.,m2=1.;
3024  LALInferenceVariableItem *item;
3025 
3026  char **info = (char **)context;
3027  char *timeID = &info[2][0];
3028 
3029  INT4 SKY_FRAME=0;
3030  if(LALInferenceCheckVariable(params,"SKY_FRAME"))
3031  SKY_FRAME=*(INT4 *)LALInferenceGetVariable(params,"SKY_FRAME");
3032 
3033  if(LALInferenceCheckVariable(params,"mass1"))
3034  {
3035  if(LALInferenceGetItem(params,"mass1")->vary != LALINFERENCE_PARAM_FIXED)
3036  {
3037  LALInferenceGetMinMaxPrior(runState->priorArgs, "mass1", (void *)&min, (void *)&max);
3038  m1 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3039  LALInferenceSetVariable(params, "mass1", &m1);
3040  }
3041  else
3042  {
3043  m1=(*(REAL8 *)LALInferenceGetVariable(params,"mass1"));
3044  }
3045  }
3046  if(LALInferenceCheckVariable(params,"mass2"))
3047  {
3048  if(LALInferenceGetItem(params,"mass2")->vary != LALINFERENCE_PARAM_FIXED)
3049  {
3050  LALInferenceGetMinMaxPrior(runState->priorArgs, "mass2", (void *)&min, (void *)&max);
3051  m2 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3052  LALInferenceSetVariable(params, "mass2", &m2);
3053  }
3054  else
3055  {
3056  m2=(*(REAL8 *)LALInferenceGetVariable(params,"mass2"));
3057  }
3058  }
3059 
3060  if(LALInferenceCheckVariable(params,"phase"))
3061  {
3062  item = LALInferenceGetItem(params,"phase");
3063  if(item->vary != LALINFERENCE_PARAM_FIXED)
3064  {
3065  LALInferenceGetMinMaxPrior(runState->priorArgs, "phase", (void *)&min, (void *)&max);
3066  double phase = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3067  LALInferenceSetVariable(params, "phase", &phase);
3068  i++;
3069  }
3070  }
3071 
3072  if(LALInferenceCheckVariable(params,"polarisation"))
3073  {
3074  item = LALInferenceGetItem(params,"polarisation");
3075  if(item->vary != LALINFERENCE_PARAM_FIXED)
3076  {
3077  LALInferenceGetMinMaxPrior(runState->priorArgs, "polarisation", (void *)&min, (void *)&max);
3078  double polarisation = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3079  LALInferenceSetVariable(params, "polarisation", &polarisation);
3080  i++;
3081  }
3082  }
3083 
3084  if(SKY_FRAME==1)
3085  {
3086  if(LALInferenceCheckVariable(params,"azimuth"))
3087  {
3088  item = LALInferenceGetItem(params,"azimuth");
3089  if(item->vary != LALINFERENCE_PARAM_FIXED)
3090  {
3091  LALInferenceGetMinMaxPrior(runState->priorArgs, "azimuth", (void *)&min, (void *)&max);
3092  double azimuth = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3093  LALInferenceSetVariable(params, "azimuth", &azimuth);
3094  i++;
3095  }
3096  }
3097  }
3098  else
3099  {
3100  if(LALInferenceCheckVariable(params,"rightascension"))
3101  {
3102  item = LALInferenceGetItem(params,"rightascension");
3103  if(item->vary != LALINFERENCE_PARAM_FIXED)
3104  {
3105  LALInferenceGetMinMaxPrior(runState->priorArgs, "rightascension", (void *)&min, (void *)&max);
3106  double rightascension = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3107  LALInferenceSetVariable(params, "rightascension", &rightascension);
3108  i++;
3109  }
3110  }
3111  }
3112 
3113  if(SKY_FRAME==1)
3114  {
3115  if(LALInferenceCheckVariable(params,"cosalpha"))
3116  {
3117  item = LALInferenceGetItem(params,"cosalpha");
3118  if(item->vary != LALINFERENCE_PARAM_FIXED)
3119  {
3120  LALInferenceGetMinMaxPrior(runState->priorArgs, "cosalpha", (void *)&min, (void *)&max);
3121  double cosalpha = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3122  LALInferenceSetVariable(params, "cosalpha", &cosalpha);
3123  i++;
3124  }
3125  }
3126  }
3127  else
3128  {
3129  if(LALInferenceCheckVariable(params,"declination"))
3130  {
3131  item = LALInferenceGetItem(params,"declination");
3132  if(item->vary != LALINFERENCE_PARAM_FIXED)
3133  {
3134  LALInferenceGetMinMaxPrior(runState->priorArgs, "declination", (void *)&min, (void *)&max);
3135  double declination = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3136  LALInferenceSetVariable(params, "declination", &declination);
3137  i++;
3138  }
3139  }
3140  }
3141 
3142  if(LALInferenceCheckVariable(params,"distance"))
3143  {
3144  item = LALInferenceGetItem(params,"distance");
3145  if(item->vary != LALINFERENCE_PARAM_FIXED)
3146  {
3147  LALInferenceGetMinMaxPrior(runState->priorArgs, "distance", (void *)&min, (void *)&max);
3148  double distance = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3149  LALInferenceSetVariable(params, "distance", &distance);
3150  i++;
3151  }
3152  }
3153 
3154  if(SKY_FRAME==1)
3155  {
3156  if(LALInferenceCheckVariable(params,"t0"))
3157  {
3158  item = LALInferenceGetItem(params,"t0");
3159  if(item->vary != LALINFERENCE_PARAM_FIXED)
3160  {
3161  LALInferenceGetMinMaxPrior(runState->priorArgs, "t0", (void *)&min, (void *)&max);
3162  double t0 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3163  LALInferenceSetVariable(params, "t0", &t0);
3164  sprintf(timeID,"%d",i);
3165  i++;
3166  }
3167  }
3168  }
3169  else
3170  {
3171  if(LALInferenceCheckVariable(params,"time"))
3172  {
3173  item = LALInferenceGetItem(params,"time");
3174  if(item->vary != LALINFERENCE_PARAM_FIXED)
3175  {
3176  LALInferenceGetMinMaxPrior(runState->priorArgs, "time", (void *)&min, (void *)&max);
3177  double tc = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3178  LALInferenceSetVariable(params, "time", &tc);
3179  sprintf(timeID,"%d",i);
3180  i++;
3181  }
3182  }
3183  }
3184 
3185  if(LALInferenceCheckVariable(params,"a_spin1"))
3186  {
3187  item = LALInferenceGetItem(params,"a_spin1");
3188  if(item->vary != LALINFERENCE_PARAM_FIXED)
3189  {
3190  LALInferenceGetMinMaxPrior(runState->priorArgs, "a_spin1", (void *)&min, (void *)&max);
3191  double a_spin1 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3192  LALInferenceSetVariable(params, "a_spin1", &a_spin1);
3193  i++;
3194  }
3195  }
3196 
3197  if(LALInferenceCheckVariable(params,"a_spin2"))
3198  {
3199  item = LALInferenceGetItem(params,"a_spin2");
3200  if(item->vary != LALINFERENCE_PARAM_FIXED)
3201  {
3202  LALInferenceGetMinMaxPrior(runState->priorArgs, "a_spin2", (void *)&min, (void *)&max);
3203  double a_spin2 = LALInferenceCubeToFlatPrior(Cube[i], min, max);
3204  LALInferenceSetVariable(params, "a_spin2", &a_spin2);
3205  i++;
3206  }
3207  }
3208 
3209  LALInferenceVariables *priorParams=runState->priorArgs;
3210  INT4 ScaleTest = LALInferenceCubeToPSDScaleParams(priorParams, params, &i, Cube, context);
3211  UINT4 ConstCalib = LALInferenceCubeToConstantCalibrationPrior(runState, params, &i, Cube, context);
3212  //UINT4 SplineCalib = LALInferenceCubeToSplineCalibrationPrior(runState, params, &i, Cube, context);
3213 
3214  if (ScaleTest==0 || ConstCalib==0 /* || SplineCalib==0 */) return 0;
3215  else return 1;
3216 }
3217 
3218 REAL8 LALInferenceFlatBoundedPrior(LALInferenceRunState *runState, LALInferenceVariables *params)
3219 {
3220  LALInferenceVariableItem *cur;
3221  REAL8 min,max;
3222  if(!params||!runState) XLAL_ERROR(XLAL_EFAULT, "Encountered NULL pointer in prior");
3223  if(!runState->priorArgs) return 0.0; /* no prior ranges specified */
3224  for(cur=params->head;cur;cur=cur->next)
3225  {
3226  if(cur->type!=LALINFERENCE_REAL8_t) continue;
3227  if(LALInferenceCheckMinMaxPrior(runState->priorArgs, cur->name))
3228  {
3229  LALInferenceGetMinMaxPrior(runState->priorArgs, cur->name, &min, &max);
3230  if (min>*(REAL8 *)cur->value || max<*(REAL8 *)cur->value) return -INFINITY;
3231  }
3232  }
3233  return 0.0;
3234 }
3235 
3236 REAL8 LALInferenceNullPrior(LALInferenceRunState UNUSED *runState, LALInferenceVariables UNUSED *params, LALInferenceModel UNUSED *model) {
3237  return 0.0;
3238 }
3239 
3240 UINT4 LALInferenceCubeToPSDScaleParams(LALInferenceVariables *priorParams, LALInferenceVariables *params, INT4 *idx, double *Cube, void UNUSED *context)
3241 {
3242  //char **info = (char **)context;
3243  //char *header = &info[1][0];
3244  //char tempstr[50];
3245 
3246  //PSD priors are Gaussian
3247  if(LALInferenceCheckVariable(params, "psdscale"))
3248  {
3249  UINT4 i;
3250  UINT4 j;
3251 
3252  REAL8 val,min,max;
3253  REAL8 mean = 1.0;
3254  UINT4 psdGaussianPrior;
3255 
3256  REAL8Vector *sigma = *((REAL8Vector **)LALInferenceGetVariable(priorParams, "psdsigma"));
3257  gsl_matrix *nparams = *((gsl_matrix **)LALInferenceGetVariable(params,"psdscale"));
3258 
3259  min=0.1;//*(REAL8 *)LALInferenceGetVariable(priorParams,"psdscale_min");
3260  max=10.0;//*(REAL8 *)LALInferenceGetVariable(priorParams,"psdscale_max");
3261 
3262  psdGaussianPrior = *(UINT4 *)LALInferenceGetVariable(priorParams,"psdGaussianPrior");
3263 
3264  for(i=0; i<(UINT4)nparams->size1; i++)
3265  {
3266  for(j=0; j<(UINT4)nparams->size2; j++)
3267  {
3268  // calculate value
3269  if (psdGaussianPrior)
3270  val = LALInferenceCubeToGaussianPrior(Cube[(*idx)],mean,sigma->data[j]);
3271  else
3272  val = LALInferenceCubeToFlatPrior(Cube[(*idx)],min,max);
3273 
3274  // set value
3275  gsl_matrix_set(nparams,i,j,val);
3276  (*idx)++;
3277  }
3278  }
3279 
3280  for(i=0; i<(UINT4)nparams->size1; i++)
3281  {
3282  for(j=0; j<(UINT4)nparams->size2; j++)
3283  {
3284  //reject prior
3285  val = gsl_matrix_get(nparams,i,j);
3286  if(val < min || val > max) return 0;
3287  }
3288  }
3289  }
3290 
3291  return 1;
3292 }
3293 
3294 /* A simple SineGaussianPrior -- will also work for other simple burst templates (gaussians) */
3295 REAL8 LALInferenceSineGaussianPrior(LALInferenceRunState *runState, LALInferenceVariables *params, LALInferenceModel *model)
3296 {
3297  REAL8 logPrior=0.0;
3298  (void)runState;
3299  LALInferenceVariableItem *item=params->head;
3300  LALInferenceVariables *priorParams=runState->priorArgs;
3301  REAL8 min, max;
3302  (void) model;
3303  /* Check boundaries */
3304  for(;item;item=item->next)
3305  {
3306  // if(item->vary!=PARAM_LINEAR || item->vary!=PARAM_CIRCULAR)
3307  if(item->vary==LALINFERENCE_PARAM_FIXED || item->vary==LALINFERENCE_PARAM_OUTPUT)
3308  continue;
3309  else
3310  {
3311  LALInferenceGetMinMaxPrior(priorParams, item->name, &min, &max);
3312  if(*(REAL8 *) item->value < min || *(REAL8 *)item->value > max) return -INFINITY;
3313  }
3314  }
3315  /*Use a distribution uniform in space volume */
3316  if(LALInferenceCheckVariable(params,"loghrss"))
3317  logPrior+=-3.0* *(REAL8 *)LALInferenceGetVariable(params,"loghrss");
3318  else if(LALInferenceCheckVariable(params,"hrss"))
3319  logPrior+=-4.0* log(*(REAL8 *)LALInferenceGetVariable(params,"hrss"));
3320  if(LALInferenceCheckVariable(params,"declination"))
3321  logPrior+=log(fabs(cos(*(REAL8 *)LALInferenceGetVariable(params,"declination"))));
3322  logPrior += LALInferenceConstantCalibrationPrior(runState, params);
3323  return(logPrior);
3324 }
3325 
3326 /**
3327  * Prior that converts from a Cube parameter in [0,1] to the flat prior bounded by x1
3328  * and x2.
3329  */
3330 REAL8 LALInferenceCubeToFlatPrior(double r, double x1, double x2)
3331 {
3332  return x1 + r * ( x2 - x1 );
3333 }
3334 
3335 /**
3336  * Prior that converts from a Cube parameter in [0,1] to the flat in log prior bounded
3337  * by x1 and x2.
3338  */
3339 REAL8 LALInferenceCubeToLogFlatPrior(double r, double x1, double x2)
3340 {
3341  double lx1, lx2;
3342  lx1 = log( x1 );
3343  lx2 = log( x2 );
3344  return exp( lx1 + r * ( lx2 - lx1 ) );
3345 }
3346 
3347 /**
3348  * Prior that converts from a Cube parameter in [0,1] to the power prior bounded by x1
3349  * and x2 with power p.
3350  */
3351 REAL8 LALInferenceCubeToPowerPrior(double p, double r, double x1, double x2)
3352 {
3353  double pp = p + 1.0;
3354  return pow(r * pow(x2, pp) + (1.0 - r) * pow(x1, pp), 1.0 / pp);
3355 }
3356 
3357 /**
3358  * Prior that converts from a Cube parameter in [0,1] to the Gaussian prior with given
3359  * mean and standard deviation.
3360  */
3361 REAL8 LALInferenceCubeToGaussianPrior(double r, double mean, double sigma)
3362 {
3363  return gsl_cdf_gaussian_Pinv(r,sigma) + mean;
3364 }
3365 
3366 /**
3367  * Prior that converts from a Cube parameter in [0,1] to the sine prior with given
3368  * min (x1) and max (x2) values
3369  */
3370 REAL8 LALInferenceCubeToSinPrior(double r, double x1, double x2)
3371 {
3372  return acos((1.0-r)*cos(x1)+r*cos(x2));
3373 }
3374 
3375 /* Functions for Fermi-Dirac prior distribution */
3376 
3377 /** \brief Return the Fermi-Dirac distribution log prior
3378  *
3379  * The function returns the log of the prior for a Fermi-Dirac distribution
3380  * \f[p(h|\sigma, r, I) = \frac{1}{\sigma\log{\left(1+e^{r} \right)}}\left(e^{((h/\sigma) - r)} + 1\right)^{-1},\f]
3381  * where \f$r = \mu/\sigma\f$ to give a more familiar form of the function. Given how it is used the function
3382  * does not actually compute the normalisation factor in the prior.
3383  */
3384 REAL8 LALInferenceFermiDiracPrior( LALInferenceVariables *priorArgs, const char *name, REAL8 value ){
3385  if ( !LALInferenceCheckFermiDiracPrior( priorArgs, name ) ){
3386  XLAL_ERROR_REAL8( XLAL_EINVAL, "No Fermi-Dirac prior given for parameter '%s'", name);
3387  }
3388 
3389  REAL8 r = 0., sigma = 0.;
3390  LALInferenceGetFermiDiracPrior(priorArgs, name, &sigma, &r);
3391 
3392  if ( value < 0. ){ return -INFINITY; } /* value must be positive */
3393  else{ return -logaddexp((value/sigma)-r, 0.); } /* log of Fermi-Dirac distribution (normalisation not required) */
3394 }
3395 
3396 /* Return the log Prior for a Gaussian Mixture Model given a value - the function returns zero
3397  * unless all required parameter values for a given GMM prior have already been passed to the
3398  * function */
3399 REAL8 LALInferenceGMMPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 value){
3400  REAL8 logPrior = 0.0;
3401 
3402  if( !LALInferenceCheckGMMPrior( priorArgs, name ) ){
3403  XLAL_ERROR_REAL8( XLAL_EINVAL, "No Gaussian Mixture Model prior given for parameter '%s'", name);
3404  }
3405 
3406  char gmmParsName[VARNAME_MAX] = "gmm_parameter_lists"; // contains list of ':'-separated lists of GMM parameters
3407 
3408  /* check if the given variable is in any of the lists */
3409  LALStringVector *parLists = *(LALStringVector **)LALInferenceGetVariable( priorArgs, gmmParsName );
3410  UINT4 numlists = parLists->length, npars = 0;
3411 
3412  UINT4 found = 0, i = 0, j = 0;
3413  TokenList *toks = NULL;
3414  for ( i = 0; i < numlists; i++ ){
3415  if ( XLALCreateTokenList( &toks, parLists->data[i], ":" ) != XLAL_SUCCESS ){ XLAL_ERROR_REAL8(XLAL_EFAILED); }
3416  for ( j = 0; j < toks->nTokens; j++ ){
3417  if ( !strcmp(toks->tokens[j], name) ){
3418  found = 1;
3419  break;
3420  }
3421  }
3422  if ( found ){ break; }
3423  }
3424 
3425  /* set value of (normalised) current parameter in priorArgs */
3426  CHAR parname[VARNAME_MAX];
3427  sprintf(parname, "%s_gmm_value", name);
3428  LALInferenceAddVariable( priorArgs, parname, &value, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED );
3429 
3430  /* check if other required parameter values have already been set */
3431  REAL8Vector *allvalues = NULL;
3432  npars = toks->nTokens;
3433  allvalues = XLALCreateREAL8Vector( npars );
3434  for ( j = 0; j < npars; j++ ){
3435  sprintf(parname, "%s_gmm_value", toks->tokens[j]);
3436  if ( LALInferenceCheckVariable( priorArgs, parname ) ){
3437  allvalues->data[j] = LALInferenceGetREAL8Variable( priorArgs, parname );
3438  }
3439  else{ /* not all values are available yet, so return zero */
3440  XLALDestroyREAL8Vector( allvalues );
3441  return 0.;
3442  }
3443  }
3444 
3445  REAL8Vector **gmmsigmas = NULL, **gmmmus = NULL, *gmmweights = NULL, *gmmdets = NULL;
3446  gsl_matrix **cor, **invcor;
3447  REAL8Vector *gmmlow = NULL, *gmmhigh = NULL;
3448  UINT4 idx = 0;
3449  CHAR *fullname = NULL;
3450 
3451  /* get GMM parameters */
3452  LALInferenceGetGMMPrior( priorArgs, name, &gmmmus, &gmmsigmas, &cor, &invcor, &gmmweights, &gmmlow, &gmmhigh, &gmmdets, &idx, &fullname );
3453 
3454  /* check values are within limits */
3455  for ( j = 0; j < npars; j++ ){
3456  if ( allvalues->data[j] < gmmlow->data[j] || allvalues->data[j] > gmmhigh->data[j] ){
3457  logPrior = -INFINITY;
3458  break;
3459  }
3460  }
3461 
3462  if ( isfinite( logPrior ) ){
3463  REAL8 thisGauss = 0.;
3464  logPrior = -INFINITY;
3465  for ( i = 0; i < gmmweights->length; i++ ){
3466  thisGauss = 0.;
3467  gsl_vector *vmu = gsl_vector_calloc( npars ); /* vector to contain value-mu */
3468  gsl_matrix *invcov = gsl_matrix_calloc( npars, npars ); /* inverse covariance matrix */
3469  gsl_matrix *tmpMat = gsl_matrix_calloc( npars, npars );
3470  for ( j = 0; j < npars; j++ ){
3471  // normalise values to be from zero mean unit variance Gaussian
3472  gsl_vector_set(vmu, j, (allvalues->data[j]-gmmmus[i]->data[j])/gmmsigmas[i]->data[j]);
3473  }
3474 
3475  /* calculate log probability */
3476  gsl_vector *tmpVec = gsl_vector_calloc( npars );
3477  gsl_blas_dgemv (CblasNoTrans, 1.0, invcor[i], vmu, 0.0, tmpVec);
3478  gsl_blas_ddot( vmu, tmpVec, &thisGauss );
3479  thisGauss *= -0.5;
3480 
3481  thisGauss += log(gmmweights->data[i]);
3482  thisGauss -= (0.5*(REAL8)npars*(LAL_LNPI + LAL_LN2) + log(gmmdets->data[i])); /* normalisation */
3483  logPrior = logaddexp(logPrior, thisGauss); /* sum Gaussianians */
3484 
3485  gsl_matrix_free( invcov );
3486  gsl_matrix_free( tmpMat );
3487  gsl_vector_free( vmu );
3488  gsl_vector_free( tmpVec );
3489  }
3490  }
3491 
3492  /* remove all values that have been set */
3493  for ( j = 0; j < npars; j++ ){
3494  sprintf(parname, "%s_gmm_value", toks->tokens[j]);
3495  LALInferenceRemoveVariable( priorArgs, parname );
3496  }
3497  XLALDestroyTokenList( toks );
3498  XLALDestroyREAL8Vector( allvalues );
3499 
3500  return logPrior;
3501 }
3502 
3503 /* Return the log Prior for a parameter that has a prior that is uniform in log space */
3504 REAL8 LALInferenceLogUniformPrior( LALInferenceVariables *priorArgs, const char *name, REAL8 value ){
3505  if ( !LALInferenceCheckLogUniformPrior( priorArgs, name ) ){
3506  XLAL_ERROR_REAL8( XLAL_EINVAL, "No log uniform prior given for parameter '%s'", name);
3507  }
3508 
3509  REAL8 min = 0., max = 0., lrat = 0.;
3510  LALInferenceGetLogUniformPrior( priorArgs, name, &min, &max );
3511 
3512  if ( value < 0. || value < min || value > max ){ return -INFINITY; }
3513  lrat = log(max/min);
3514 
3515  return -log(value*lrat);
3516 }
XLALCheckBurstApproximantFromString
int XLALCheckBurstApproximantFromString(const CHAR *inString)
Definition: LALInferenceBurstRoutines.c:101
mean
static REAL8 mean(REAL8 *array, int N)
Definition: LALInferenceNestedSampler.c:387
LALINFERENCE_PRIOR_MIN
#define LALINFERENCE_PRIOR_MIN(x, y)
Definition: LALInferencePrior.c:1853
outerIntegrand
static double outerIntegrand(double M1, void *voData)
Definition: LALInferencePrior.c:1855
qInnerIntegrand
static double qInnerIntegrand(double M2, void *viData)
Definition: LALInferencePrior.c:1813
LALInferenceCubeToConstantCalibrationPrior
UINT4 LALInferenceCubeToConstantCalibrationPrior(LALInferenceRunState *runState, LALInferenceVariables *params, INT4 *idx, double *Cube, void UNUSED *context)
Definition: LALInferencePrior.c:250
LALInferenceInspiralSkyLocPrior
REAL8 LALInferenceInspiralSkyLocPrior(LALInferenceRunState *runState, LALInferenceVariables *params, UNUSED LALInferenceModel *model)
Definition: LALInferencePrior.c:1247
LALInferenceAnalyticNullPrior
REAL8 LALInferenceAnalyticNullPrior(LALInferenceRunState UNUSED *runState, LALInferenceVariables *params, LALInferenceModel UNUSED *model)
Definition: LALInferencePrior.c:2985
LALInferenceGlitchPrior
static REAL8 LALInferenceGlitchPrior(LALInferenceRunState *runState, LALInferenceVariables *params)
Definition: LALInferencePrior.c:319
etaInnerIntegrand
static double etaInnerIntegrand(double M2, void *viData)
Definition: LALInferencePrior.c:1826
REAL8max
static REAL8 REAL8max(REAL8 a, REAL8 b)
Definition: LALInferencePrior.c:47
LALInferenceCubeToPSDScaleParams
UINT4 LALInferenceCubeToPSDScaleParams(LALInferenceVariables *priorParams, LALInferenceVariables *params, INT4 *idx, double *Cube, void UNUSED *context)
Definition: LALInferencePrior.c:3240
LALInferenceConstantCalibrationPrior
static REAL8 LALInferenceConstantCalibrationPrior(LALInferenceRunState *runState, LALInferenceVariables *params)
Definition: LALInferencePrior.c:201
LALInferenceAnalyticCubeToPrior
UINT4 LALInferenceAnalyticCubeToPrior(LALInferenceRunState *runState, LALInferenceVariables *params, UNUSED LALInferenceModel *model, double *Cube, void *context)
Definition: LALInferencePrior.c:3020
LALInferenceNullPrior
REAL8 LALInferenceNullPrior(LALInferenceRunState UNUSED *runState, LALInferenceVariables UNUSED *params, LALInferenceModel UNUSED *model)
Definition: LALInferencePrior.c:3236
LALINFERENCE_PRIOR_SQR
#define LALINFERENCE_PRIOR_SQR(x)
Definition: LALInferencePrior.c:1824
LALInferenceInspiralSkyLocCubeToPrior
UINT4 LALInferenceInspiralSkyLocCubeToPrior(LALInferenceRunState *runState, LALInferenceVariables *params, UNUSED LALInferenceModel *model, double *Cube, void *context)
Definition: LALInferencePrior.c:1408
LALInferencePSDPrior
static REAL8 LALInferencePSDPrior(LALInferenceRunState *runState, LALInferenceVariables *params)
Definition: LALInferencePrior.c:391
max
#define max(a, b)
Definition: LALInferenceRemoveLines.c:80
ppt
ProcessParamsTable * ppt
Definition: LALInferenceTest.c:379
j
int j
Definition: LALInferenceTest.c:381
ptr
ProcessParamsTable * ptr
Definition: LALInferenceTest.c:379
MMax
#define MMax
c1
const double c1
c2
const double c2
c0
const double c0
delta
static double double delta
fprintf
#define fprintf
XLALDestroyTokenList
void XLALDestroyTokenList(TokenList *list)
XLALCreateTokenList
int XLALCreateTokenList(TokenList **list, const CHAR *string, const CHAR *delimiters)
name
const char *const name
XLAL_CALLGSL
#define XLAL_CALLGSL(statement)
s
int s
e
double e
pp
const double pp
Q
const double Q
data
sigmaKerr data[0]
LAL_2_SQRTPI
#define LAL_2_SQRTPI
LAL_LN2
#define LAL_LN2
LAL_LNPI
#define LAL_LNPI
LAL_PI
#define LAL_PI
LAL_TWOPI
#define LAL_TWOPI
LAL_SQRT1_2
#define LAL_SQRT1_2
UINT8
uint64_t UINT8
REAL8
double REAL8
INT8
int64_t INT8
CHAR
char CHAR
UINT4
uint32_t UINT4
INT4
int32_t INT4
REAL4
float REAL4
LALInferenceMcQ2Masses
void LALInferenceMcQ2Masses(double mc, double q, double *m1, double *m2)
Convert from Mc, q space to m1, m2 space (q = m2/m1, with m1 > m2).
Definition: LALInference.c:2304
VARNAME_MAX
#define VARNAME_MAX
Definition: LALInference.h:50
LALInferenceGetUINT4Variable
UINT4 LALInferenceGetUINT4Variable(LALInferenceVariables *vars, const char *name)
Definition: LALInference.c:3993
LALInferenceAddVariable
void LALInferenceAddVariable(LALInferenceVariables *vars, const char *name, const void *value, LALInferenceVariableType type, LALInferenceParamVaryType vary)
Add a variable named name to vars with initial value referenced by value.
Definition: LALInference.c:395
LALInferenceGetREAL8Variable
REAL8 LALInferenceGetREAL8Variable(LALInferenceVariables *vars, const char *name)
Definition: LALInference.c:4039
LALInferenceEOSPhysicalCheck
int LALInferenceEOSPhysicalCheck(LALInferenceVariables *params, ProcessParamsTable *commandLine)
Check for causality violation and mass conflict given masses and eos.
Definition: LALInference.c:2410
LALInferenceRemoveVariable
void LALInferenceRemoveVariable(LALInferenceVariables *vars, const char *name)
Remove name from vars Frees the memory for the name structure and its contents.
Definition: LALInference.c:450
LALInferenceParamVaryType
LALInferenceParamVaryType
An enumerated type for denoting the topology of a parameter.
Definition: LALInference.h:127
LALInferenceGetINT4Variable
INT4 LALInferenceGetINT4Variable(LALInferenceVariables *vars, const char *name)
Definition: LALInference.c:3947
LALInferenceKDDrawEigenFrame
void LALInferenceKDDrawEigenFrame(gsl_rng *rng, LALInferenceKDTree *tree, REAL8 *pt, size_t Npts)
Draws a pt uniformly from a randomly chosen cell of tree.
Definition: LALInference.c:3167
LALInferenceCheckPositiveDefinite
UINT4 LALInferenceCheckPositiveDefinite(gsl_matrix *matrix, UINT4 dim)
Check matrix is positive definite.
Definition: LALInference.c:3281
LALInferenceMcEta2Masses
void LALInferenceMcEta2Masses(double mc, double eta, double *m1, double *m2)
Convert from Mc, eta space to m1, m2 space (note m1 > m2).
Definition: LALInference.c:2292
LALInferenceGetProcParamVal
ProcessParamsTable * LALInferenceGetProcParamVal(ProcessParamsTable *procparams, const char *name)
Returns the element of the process params table with "name".
Definition: LALInference.c:1762
LALInferenceGetVariable
void * LALInferenceGetVariable(const LALInferenceVariables *vars, const char *name)
Return a pointer to the memory the variable vars is stored in specified by name User must cast this p...
Definition: LALInference.c:238
XLALMultiNormalDeviates
void XLALMultiNormalDeviates(REAL4Vector *vector, gsl_matrix *matrix, UINT4 dim, RandomParams *randParam)
Draw a random multivariate vector from Gaussian distr given covariance matrix.
Definition: LALInference.c:3323
LALInferenceCheckVariable
int LALInferenceCheckVariable(LALInferenceVariables *vars, const char *name)
Checks for name being present in vars returns 1(==true) or 0.
Definition: LALInference.c:525
LALInferenceVariableType
LALInferenceVariableType
An enumerated type for denoting the type of a variable.
Definition: LALInference.h:104
LALInferenceGetVariableVaryType
LALInferenceParamVaryType LALInferenceGetVariableVaryType(LALInferenceVariables *vars, const char *name)
Get the LALInferenceParamVaryType of the parameter named name in vars see the declaration of LALInfer...
Definition: LALInference.c:226
LALInferenceGetItem
LALInferenceVariableItem * LALInferenceGetItem(const LALInferenceVariables *vars, const char *name)
Return the list node for "name" - do not rely on this.
Definition: LALInference.c:175
LALInferenceSetVariable
void LALInferenceSetVariable(LALInferenceVariables *vars, const char *name, const void *value)
Set a variable named name in vars with a value.
Definition: LALInference.c:352
LALInferenceCheckVariableNonFixed
int LALInferenceCheckVariableNonFixed(LALInferenceVariables *vars, const char *name)
Checks for name being present in vars and having type LINEAR or CIRCULAR.
Definition: LALInference.c:503
LALInferenceKDREAL8ToVariables
void LALInferenceKDREAL8ToVariables(LALInferenceVariables *params, REAL8 *pt, LALInferenceVariables *templt)
Fills in the non-fixed variables in params from the given REAL8 array.
Definition: LALInference.c:3155
LALINFERENCE_PARAM_OUTPUT
@ LALINFERENCE_PARAM_OUTPUT
A parameter that never changes, functions should respect this.
Definition: LALInference.h:131
LALINFERENCE_PARAM_FIXED
@ LALINFERENCE_PARAM_FIXED
A parameter that is cyclic, such as an angle between 0 and 2pi.
Definition: LALInference.h:130
LALINFERENCE_PARAM_CIRCULAR
@ LALINFERENCE_PARAM_CIRCULAR
A parameter that simply has a maximum and a minimum.
Definition: LALInference.h:129
LALINFERENCE_PARAM_LINEAR
@ LALINFERENCE_PARAM_LINEAR
Definition: LALInference.h:128
LALINFERENCE_REAL8_t
@ LALINFERENCE_REAL8_t
Definition: LALInference.h:109
LALINFERENCE_REAL4_t
@ LALINFERENCE_REAL4_t
Definition: LALInference.h:108
LALINFERENCE_INT4_t
@ LALINFERENCE_INT4_t
Definition: LALInference.h:105
LALINFERENCE_REAL8Vector_t
@ LALINFERENCE_REAL8Vector_t
Definition: LALInference.h:113
LALINFERENCE_void_ptr_t
@ LALINFERENCE_void_ptr_t
Definition: LALInference.h:119
LALINFERENCE_UINT4_t
@ LALINFERENCE_UINT4_t
Definition: LALInference.h:107
LALINFERENCE_gslMatrix_t
@ LALINFERENCE_gslMatrix_t
Definition: LALInference.h:112
LALINFERENCE_INT8_t
@ LALINFERENCE_INT8_t
Definition: LALInference.h:106
LALInferenceNetworkSNR
void LALInferenceNetworkSNR(LALInferenceVariables *currentParams, LALInferenceIFOData *data, LALInferenceModel *model)
Calculate the SNR across the network.
Definition: LALInferenceLikelihood.c:1866
LALInferenceRemoveCorrelatedPrior
void LALInferenceRemoveCorrelatedPrior(LALInferenceVariables *priorArgs)
Remove the correlation coefficient matrix and index for a parameter from the priorArgs list.
Definition: LALInferencePrior.c:2146
LALInferenceGMMPrior
REAL8 LALInferenceGMMPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 value)
Calculate the log probability for the Gaussian Mixture Model prior.
Definition: LALInferencePrior.c:3399
LALInferenceAddGaussianPrior
void LALInferenceAddGaussianPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *mu, REAL8 *sigma, LALInferenceVariableType type)
Function to add the mu and sigma values for the Gaussian prior onto the priorArgs.
Definition: LALInferencePrior.c:2003
LALInferenceComputePriorMassNorm
REAL8 LALInferenceComputePriorMassNorm(const double MMin, const double MMax, const double MTotMax, const double McMin, const double McMax, const double massRatioMin, const double massRatioMax, const char *massRatioName)
Computes the numerical normalization of the mass prior applying all cuts in the mass plane implied b...
Definition: LALInferencePrior.c:1880
LALInferenceInspiralPrior
REAL8 LALInferenceInspiralPrior(LALInferenceRunState *runState, LALInferenceVariables *params, LALInferenceModel *model)
Return the logarithmic prior density of the variables specified, for the non-spinning/spinning inspir...
Definition: LALInferencePrior.c:442
LALInferenceFlatBoundedPrior
REAL8 LALInferenceFlatBoundedPrior(LALInferenceRunState *runState, LALInferenceVariables *params)
Prior that checks for minimum and maximum prior range specified in runState->priorArgs and returns 0....
Definition: LALInferencePrior.c:3218
LALInferenceCheckMinMaxPrior
int LALInferenceCheckMinMaxPrior(LALInferenceVariables *priorArgs, const char *name)
Check for types of standard prior.
Definition: LALInferencePrior.c:1963
LALInferenceCubeToPowerPrior
REAL8 LALInferenceCubeToPowerPrior(double p, double r, double x1, double x2)
Prior that converts from a Cube parameter in [0,1] to the power prior bounded by x1 and x2 with power...
Definition: LALInferencePrior.c:3351
LALInferenceAddGMMPrior
void LALInferenceAddGMMPrior(LALInferenceVariables *priorArgs, const char *name, REAL8Vector ***mus, gsl_matrix ***covs, REAL8Vector **weights, REAL8Vector **minrange, REAL8Vector **maxrange)
Add a Gaussian Mixture Model prior.
Definition: LALInferencePrior.c:2203
LALInferenceDrawNameFromPrior
void LALInferenceDrawNameFromPrior(LALInferenceVariables *output, LALInferenceVariables *priorArgs, char *name, LALInferenceVariableType type, gsl_rng *rdm)
Draw an individual variable from its prior range.
Definition: LALInferencePrior.c:2736
LALInferenceAddMinMaxPrior
void LALInferenceAddMinMaxPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *min, REAL8 *max, LALInferenceVariableType type)
Function to add the minimum and maximum values for the uniform prior onto the priorArgs.
Definition: LALInferencePrior.c:1934
LALInferenceInitCBCPrior
void LALInferenceInitCBCPrior(LALInferenceRunState *runState)
Initialize the prior based on command line arguments.
Definition: LALInferencePrior.c:52
LALInferenceRemoveFermiDiracPrior
void LALInferenceRemoveFermiDiracPrior(LALInferenceVariables *priorArgs, const char *name)
Function to remove the r and sigma values for the Fermi-Dirac prior onto the priorArgs.
Definition: LALInferencePrior.c:2576
LALInferenceRemoveMinMaxPrior
void LALInferenceRemoveMinMaxPrior(LALInferenceVariables *priorArgs, const char *name)
Function to remove the minimum and maximum values for the uniform prior onto the priorArgs.
Definition: LALInferencePrior.c:1950
LALInferenceCubeToLogFlatPrior
REAL8 LALInferenceCubeToLogFlatPrior(double r, double x1, double x2)
Prior that converts from a Cube parameter in [0,1] to the flat in log prior bounded by x1 and x2.
Definition: LALInferencePrior.c:3339
LALInferenceSineGaussianPrior
REAL8 LALInferenceSineGaussianPrior(LALInferenceRunState *runState, LALInferenceVariables *params, LALInferenceModel *model)
Definition: LALInferencePrior.c:3295
LALInferenceCubeToSinPrior
REAL8 LALInferenceCubeToSinPrior(double r, double x1, double x2)
Prior that converts from a Cube parameter in [0,1] to the sine prior with given min (x1) and max (x2)...
Definition: LALInferencePrior.c:3370
LALInferenceAddLogUniformPrior
void LALInferenceAddLogUniformPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *xmin, REAL8 *xmax, LALInferenceVariableType type)
Add a log-uniform prior.
Definition: LALInferencePrior.c:2626
LALInferenceRemoveGaussianPrior
void LALInferenceRemoveGaussianPrior(LALInferenceVariables *priorArgs, const char *name)
Function to remove the mu and sigma values for the Gaussian prior onto the priorArgs.
Definition: LALInferencePrior.c:2018
LALInferenceLogUniformPrior
REAL8 LALInferenceLogUniformPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 value)
Definition: LALInferencePrior.c:3504
LALInferenceCyclicReflectiveBound
void LALInferenceCyclicReflectiveBound(LALInferenceVariables *parameter, LALInferenceVariables *priorArgs)
Apply cyclic and reflective boundaries to parameter to bring it back within the allowed prior ranges ...
Definition: LALInferencePrior.c:1113
LALInferenceRotateInitialPhase
void LALInferenceRotateInitialPhase(LALInferenceVariables *parameter)
Rotate initial phase if polarisation angle is cyclic around ranges.
Definition: LALInferencePrior.c:1204
LALInferenceCheckCorrelatedPrior
int LALInferenceCheckCorrelatedPrior(LALInferenceVariables *priorArgs, const char *name)
Check for the existance of a correlation coefficient matrix and index for a parameter from the priorA...
Definition: LALInferencePrior.c:2176
LALInferenceDrawFromPrior
void LALInferenceDrawFromPrior(LALInferenceVariables *output, LALInferenceVariables *priorArgs, gsl_rng *rdm)
Draw variables from the prior ranges.
Definition: LALInferencePrior.c:2690
LALInferenceInitLIBPrior
void LALInferenceInitLIBPrior(LALInferenceRunState *runState)
Initialize the LIB prior based on command line arguments.
Definition: LALInferencePrior.c:177
LALInferenceFermiDiracPrior
REAL8 LALInferenceFermiDiracPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 value)
Return the Fermi-Dirac distribution log prior.
Definition: LALInferencePrior.c:3384
LALInferenceCheckGMMPrior
int LALInferenceCheckGMMPrior(LALInferenceVariables *priorArgs, const char *name)
Check for a Gaussian Mixture Model prior.
Definition: LALInferencePrior.c:2424
LALInferenceGetLogUniformPrior
void LALInferenceGetLogUniformPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *xmin, REAL8 *xmax)
Get the xmin and xmax values of the log-uniform prior from the priorArgs list, given a name.
Definition: LALInferencePrior.c:2663
LALInferenceGetMinMaxPrior
void LALInferenceGetMinMaxPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *min, REAL8 *max)
Get the minimum and maximum values of the uniform prior from the priorArgs list, given a name.
Definition: LALInferencePrior.c:1974
LALInferenceCheckGaussianPrior
int LALInferenceCheckGaussianPrior(LALInferenceVariables *priorArgs, const char *name)
Check for a Gaussian prior (with a mean and variance)
Definition: LALInferencePrior.c:1993
LALInferenceGetGMMPrior
void LALInferenceGetGMMPrior(LALInferenceVariables *priorArgs, const char *name, REAL8Vector ***mus, REAL8Vector ***sigmas, gsl_matrix ***cors, gsl_matrix ***invcors, REAL8Vector **weights, REAL8Vector **minrange, REAL8Vector **maxrange, REAL8Vector **dets, UINT4 *idx, CHAR **fullname)
Get the parameters defining a Gaussian Mixture Model prior.
Definition: LALInferencePrior.c:2341
LALInferenceCheckFermiDiracPrior
int LALInferenceCheckFermiDiracPrior(LALInferenceVariables *priorArgs, const char *name)
Check for a Fermi-Dirac prior (with a r and sigma parameter)
Definition: LALInferencePrior.c:2551
LALInferenceCheckLogUniformPrior
int LALInferenceCheckLogUniformPrior(LALInferenceVariables *priorArgs, const char *name)
Check for a log-uniform prior (with xmin and xmax parameters)
Definition: LALInferencePrior.c:2612
within_malmquist
UINT4 within_malmquist(LALInferenceRunState *runState, LALInferenceVariables *params, LALInferenceModel *model)
Definition: LALInferencePrior.c:2950
LALInferenceGetGaussianPrior
void LALInferenceGetGaussianPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *mu, REAL8 *sigma)
Get the mu and sigma values of the Gaussian prior from the priorArgs list, given a name.
Definition: LALInferencePrior.c:2032
LALInferenceCubeToGaussianPrior
REAL8 LALInferenceCubeToGaussianPrior(double r, double mean, double sigma)
Prior that converts from a Cube parameter in [0,1] to the Gaussian prior with given mean and standard...
Definition: LALInferencePrior.c:3361
logGlitchAmplitudeDensity
REAL8 logGlitchAmplitudeDensity(REAL8 A, REAL8 Q, REAL8 f)
Return the log Prior for the glitch amplitude.
Definition: LALInferencePrior.c:305
LALInferenceRemoveGMMPrior
void LALInferenceRemoveGMMPrior(LALInferenceVariables *priorArgs, const char *name)
Remove a Gaussian Mixture Model prior.
Definition: LALInferencePrior.c:2485
LALInferenceInspiralCubeToPrior
UINT4 LALInferenceInspiralCubeToPrior(LALInferenceRunState *runState, LALInferenceVariables *params, LALInferenceModel *model, double *Cube, void *context)
Convert the hypercube parameter to physical parameters, for the non-spinning/spinning inspiral signal...
Definition: LALInferencePrior.c:707
LALInferenceAddCorrelatedPrior
void LALInferenceAddCorrelatedPrior(LALInferenceVariables *priorArgs, const char *name, gsl_matrix **cor, REAL8 *mu, REAL8 *sigma, UINT4 *idx)
Function to add a correlation matrix and parameter index for a prior defined as part of a multivariat...
Definition: LALInferencePrior.c:2054
LALInferenceRemoveLogUniformPrior
void LALInferenceRemoveLogUniformPrior(LALInferenceVariables *priorArgs, const char *name)
Function to remove the min and max values for the log-uniform prior from the priorArgs.
Definition: LALInferencePrior.c:2648
LALInferenceCubeToFlatPrior
REAL8 LALInferenceCubeToFlatPrior(double r, double x1, double x2)
Prior that converts from a Cube parameter in [0,1] to the flat prior bounded by x1 and x2.
Definition: LALInferencePrior.c:3330
LALInferenceAddFermiDiracPrior
void LALInferenceAddFermiDiracPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *sigma, REAL8 *r, LALInferenceVariableType type)
Add a Fermi-Dirac prior.
Definition: LALInferencePrior.c:2561
LALInferenceGetCorrelatedPrior
void LALInferenceGetCorrelatedPrior(LALInferenceVariables *priorArgs, const char *name, gsl_matrix **cor, gsl_matrix **invcor, REAL8 *mu, REAL8 *sigma, UINT4 *idx)
Get the correlation coefficient matrix and index for a parameter from the priorArgs list.
Definition: LALInferencePrior.c:2106
LALInferenceGetFermiDiracPrior
void LALInferenceGetFermiDiracPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *sigma, REAL8 *r)
Get the r and sigma values of the Fermi-Dirac prior from the priorArgs list, given a name.
Definition: LALInferencePrior.c:2589
XLALCalloc
void * XLALCalloc(size_t m, size_t n)
r
static const INT4 r
m
static const INT4 m
q
static const INT4 q
XLALCreateRandomParams
RandomParams * XLALCreateRandomParams(INT4 seed)
XLALDestroyRandomParams
void XLALDestroyRandomParams(RandomParams *params)
a
static const INT4 a
XLALAppendString2Vector
LALStringVector * XLALAppendString2Vector(LALStringVector *vect, const CHAR *string)
XLALDestroyREAL4Vector
void XLALDestroyREAL4Vector(REAL4Vector *vector)
XLALCreateREAL4Vector
REAL4Vector * XLALCreateREAL4Vector(UINT4 length)
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLAL_ERROR_VOID
#define XLAL_ERROR_VOID(...)
XLAL_ERROR_REAL8
#define XLAL_ERROR_REAL8(...)
XLAL_ERROR
#define XLAL_ERROR(...)
XLALPrintWarning
int int XLALPrintWarning(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_EDATA
XLAL_EDATA
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFAULT
XLAL_ENAME
XLAL_ENAME
XLAL_EFUNC
XLAL_EFUNC
XLAL_EDOM
XLAL_EDOM
XLAL_EINVAL
XLAL_EINVAL
XLAL_EFAILED
XLAL_EFAILED
XLAL_FAILURE
XLAL_FAILURE
logaddexp
static double logaddexp(double x, double y)
Definition: logaddexp.h:37
mu
list mu
cbcBayesBurstPostProc.help
help
Definition: cbcBayesBurstPostProc.py:1198
cbcBayesBurstPostProc.type
type
Definition: cbcBayesBurstPostProc.py:1204
lalinference.git_version.status
string status
Definition: git_version.py:32
lalinference.tiger.make_injtimes.dist
dist
Definition: make_injtimes.py:581
lalinference.tiger.make_injtimes.n
n
Definition: make_injtimes.py:576
lalinference_compute_roq_weights.i
int i
Definition: lalinference_compute_roq_weights.py:269
lalinference_compute_roq_weights.tc
tc
Definition: lalinference_compute_roq_weights.py:335
lalinference_evolve_spins_and_append_samples.phi12
phi12
Definition: lalinference_evolve_spins_and_append_samples.py:113
lalinference_merge_posteriors.weights
list weights
Definition: lalinference_merge_posteriors.py:40
lalinference_multi_pipe.cp
cp
Definition: lalinference_multi_pipe.py:82
lalinference_review_test.output
output
Definition: lalinference_review_test.py:76
t0
double t0
LALInferenceIFOData
Structure to contain IFO data.
Definition: LALInference.h:625
LALInferenceIFOData::name
char name[DETNAMELEN]
Definition: LALInference.h:626
LALInferenceIFOData::next
struct tagLALInferenceIFOData * next
ROQ data.
Definition: LALInference.h:659
LALInferenceKDTree
The kD trees in LALInference are composed of cells.
Definition: LALInference.h:915
LALInferenceKDTree::dim
size_t dim
Size of the pts buffer.
Definition: LALInference.h:918
LALInferenceModel
Structure to constain a model and its parameters.
Definition: LALInference.h:436
LALInferenceModel::ifo_SNRs
REAL8 * ifo_SNRs
Array of single-IFO likelihoods at params
Definition: LALInference.h:446
LALInferenceModel::SNR
REAL8 SNR
Likelihood value at params
Definition: LALInference.h:444
LALInferenceRunState
Structure containing inference run state This includes pointers to the function types required to run...
Definition: LALInference.h:592
LALInferenceRunState::commandLine
ProcessParamsTable * commandLine
Definition: LALInference.h:593
LALInferenceRunState::CubeToPrior
LALInferenceCubeToPriorFunction CubeToPrior
The prior for the parameters.
Definition: LALInference.h:598
LALInferenceRunState::data
struct tagLALInferenceIFOData * data
Log sample, i.e.
Definition: LALInference.h:601
LALInferenceRunState::priorArgs
LALInferenceVariables * priorArgs
Common arguments needed by proposals, to be copied to thread->cycle.
Definition: LALInference.h:603
LALInferenceRunState::prior
LALInferencePriorFunction prior
The algorithm's single iteration function.
Definition: LALInference.h:597
LALInferenceVariableItem
The LALInferenceVariableItem list node structure This should only be accessed using the accessor func...
Definition: LALInference.h:154
LALInferenceVariableItem::name
char name[VARNAME_MAX]
Definition: LALInference.h:155
LALInferenceVariableItem::type
LALInferenceVariableType type
Definition: LALInference.h:157
LALInferenceVariableItem::next
struct tagVariableItem * next
Definition: LALInference.h:159
LALInferenceVariableItem::vary
LALInferenceParamVaryType vary
Definition: LALInference.h:158
LALInferenceVariables
The LALInferenceVariables structure to contain a set of parameters Implemented as a linked list of LA...
Definition: LALInference.h:170
LALInferenceVariables::head
LALInferenceVariableItem * head
Definition: LALInference.h:171
ProcessParamsTable::value
CHAR value[LIGOMETA_VALUE_MAX]
REAL4Vector::data
REAL4 * data
REAL8Vector::data
REAL8 * data
TokenList::tokens
CHAR ** tokens
TokenList::nTokens
UINT4 nTokens
UINT4Vector::data
UINT4 * data
outerData::innerIntegrand
gsl_function innerIntegrand
Definition: LALInferencePrior.c:1850
outerData::wsInner
gsl_integration_workspace * wsInner
Definition: LALInferencePrior.c:1840
V
double V