LALInferenceLikelihood.c 79.3 KB
Newer Older
1
/*
2 3
 *  LALInferenceLikelihood.c:  Bayesian Followup likelihood functions
 *
Will Farr's avatar
Will Farr committed
4 5
 *  Copyright (C) 2009 Ilya Mandel, Vivien Raymond, Christian Roever,
 *  Marc van der Sluys and John Veitch, Will M. Farr
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
 *
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with with program; see the file COPYING. If not, write to the
 *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
 *  MA  02111-1307  USA
 */

24
#include <complex.h>
25
#include <assert.h>
26
#include <lal/LALInferenceLikelihood.h>
27
#include <lal/LALInferencePrior.h>
28
#include <lal/LALInference.h>
29 30 31 32 33 34 35
#include <lal/DetResponse.h>
#include <lal/TimeDelay.h>
#include <lal/TimeSeries.h>
#include <lal/Units.h>
#include <lal/Sequence.h>
#include <lal/FrequencySeries.h>
#include <lal/TimeFreqFFT.h>
36
#include <lal/LALInferenceDistanceMarg.h>
37

38
#include <gsl/gsl_sf_bessel.h>
39 40
#include <gsl/gsl_sf_dawson.h>
#include <gsl/gsl_sf_erf.h>
41
#include <gsl/gsl_complex_math.h>
42
#include <lal/LALInferenceTemplate.h>
43

44 45
#include "logaddexp.h"

46
#include <lal/distance_integrator.h>
47

48 49
typedef enum
{
50 51 52 53 54 55
  GAUSSIAN = 1,
  STUDENTT = 2,
  MARGPHI = 4,
  MARGTIME = 8,
  MARGTIMEPHI = 16,
  MARGDIST = 32
56 57
} LALInferenceLikelihoodFlags;

58

59 60 61 62 63 64
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

John Douglas Veitch's avatar
John Douglas Veitch committed
65 66 67 68
#ifndef _OPENMP
#define omp ignore
#endif

69 70 71 72 73
static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *currentParams,
                                               LALInferenceIFOData *data,
                                               LALInferenceModel *model,
                                               LALInferenceLikelihoodFlags marginalisationflags);

74
static double integrate_interpolated_log(double h, REAL8 *log_ys, size_t n, double *imean, size_t *imax);
75

76 77 78 79 80 81 82 83 84 85 86 87 88
static int get_calib_spline(LALInferenceVariables *vars, const char *ifoname, REAL8Vector **logfreqs, REAL8Vector **amps, REAL8Vector **phases);
static int get_calib_spline(LALInferenceVariables *vars, const char *ifoname, REAL8Vector **logfreqs, REAL8Vector **amps, REAL8Vector **phases)
{
  UINT4 npts = LALInferenceGetUINT4Variable(vars, "spcal_npts");
  char ampname[VARNAME_MAX];
  char phasename[VARNAME_MAX];
  char freqname[VARNAME_MAX];
  if(!*logfreqs) *logfreqs = XLALCreateREAL8Vector(npts);
  if(!*amps) *amps = XLALCreateREAL8Vector(npts);
  if(!*phases) *phases = XLALCreateREAL8Vector(npts);
  assert((*logfreqs)->length==npts);
  assert((*amps)->length==npts);
  assert((*phases)->length==npts);
89

90 91
  for(UINT4 i=0;i<npts;i++)
  {
92 93 94
    if((VARNAME_MAX <= snprintf(freqname, VARNAME_MAX, "%s_spcal_logfreq_%i", ifoname, i))) XLAL_ERROR(XLAL_EINVAL,"Variable name too long");
    if((VARNAME_MAX <= snprintf(ampname, VARNAME_MAX, "%s_spcal_amp_%i", ifoname, i))) XLAL_ERROR(XLAL_EINVAL,"Variable name too long");
    if((VARNAME_MAX <= snprintf(phasename, VARNAME_MAX, "%s_spcal_phase_%i", ifoname, i))) XLAL_ERROR(XLAL_EINVAL,"Variable name too long");
95 96 97 98 99 100 101

    (*logfreqs)->data[i] = LALInferenceGetREAL8Variable(vars, freqname);
    (*amps)->data[i] =  LALInferenceGetREAL8Variable(vars, ampname);
    (*phases)->data[i] = LALInferenceGetREAL8Variable(vars, phasename);
  }
  return(XLAL_SUCCESS);
}
102

103 104 105
void LALInferenceInitLikelihood(LALInferenceRunState *runState)
{
    char help[]="\
Ben Farr's avatar
Ben Farr committed
106 107 108 109 110
    ----------------------------------------------\n\
    --- Likelihood Arguments ---------------------\n\
    ----------------------------------------------\n\
    (--zeroLogLike)                  Use flat, null likelihood\n\
    (--studentTLikelihood)           Use the Student-T Likelihood that marginalizes over noise\n\
111 112
    (--correlatedGaussianLikelihood) Use analytic, correlated Gaussian for Likelihood Z=-21.3\n\
    (--bimodalGaussianLikelihood)    Use analytic, bimodal correlated Gaussian for Likelihood Z=-25.9\n\
Ben Farr's avatar
Ben Farr committed
113 114 115 116 117
    (--rosenbrockLikelihood)         Use analytic, Rosenbrock banana for Likelihood\n\
    (--noiseonly)                    Using noise-only likelihood\n\
    (--margphi)                      Using marginalised phase likelihood\n\
    (--margtime)                     Using marginalised time likelihood\n\
    (--margtimephi)                  Using marginalised in time and phase likelihood\n\
John Douglas Veitch's avatar
John Douglas Veitch committed
118 119
    (--margdist)                     Using marginalisation in distance with d^2 prior (compatible with --margphi and --margtimephi)\n\
    (--margdist-comoving)            Using marginalisation in distance with uniform-in-comoving-volume prior (compatible with --margphi and --margtimephi)\n\
Ben Farr's avatar
Ben Farr committed
120
    \n";
121 122

    /* Print command line arguments if help requested */
Vivien Raymond's avatar
Vivien Raymond committed
123 124
    LALInferenceIFOData *ifo=NULL;
    if(runState == NULL || LALInferenceGetProcParamVal(runState->commandLine,"--help"))
125 126
    {
        fprintf(stdout,"%s",help);
Vivien Raymond's avatar
Vivien Raymond committed
127 128 129 130 131 132
        if (runState){
            ifo=runState->data;
            while(ifo) {
                fprintf(stdout,"(--dof-%s DoF)\tDegrees of freedom for %s\n",ifo->name,ifo->name);
                ifo=ifo->next;
            }
133 134 135 136
        }
        return;
    }

Vivien Raymond's avatar
Vivien Raymond committed
137 138 139 140 141 142 143
    ProcessParamsTable *commandLine=runState->commandLine;
    ifo=runState->data;

    LALInferenceThreadState *thread = runState->threads[0];

    REAL8 nullLikelihood = 0.0; // Populated if such a thing exists

144 145 146 147 148 149 150 151 152 153 154 155 156 157
   if (LALInferenceGetProcParamVal(commandLine, "--zeroLogLike")) {
    /* Use zero log(L) */
    runState->likelihood=&LALInferenceZeroLogLikelihood;
   } else if (LALInferenceGetProcParamVal(commandLine, "--correlatedGaussianLikelihood")) {
    runState->likelihood=&LALInferenceCorrelatedAnalyticLogLikelihood;
   } else if (LALInferenceGetProcParamVal(commandLine, "--bimodalGaussianLikelihood")) {
    runState->likelihood=&LALInferenceBimodalCorrelatedAnalyticLogLikelihood;
   } else if (LALInferenceGetProcParamVal(commandLine, "--rosenbrockLikelihood")) {
    runState->likelihood=&LALInferenceRosenbrockLogLikelihood;
   } else if (LALInferenceGetProcParamVal(commandLine, "--studentTLikelihood")) {
    fprintf(stderr, "Using Student's T Likelihood.\n");
    runState->likelihood=&LALInferenceFreqDomainStudentTLogLikelihood;

    /* Set the noise model evidence to the student t model value */
158 159 160 161 162 163 164
    LALInferenceTemplateNullFreqdomain(thread->model);
    LALInferenceTemplateFunction temp = thread->model->templt;
    thread->model->templt = &LALInferenceTemplateNullFreqdomain;
    REAL8 noiseZ = LALInferenceFreqDomainStudentTLogLikelihood(thread->currentParams, runState->data, thread->model);
    thread->model->templt = temp;
    LALInferenceAddVariable(runState->algorithmParams, "logZnoise", &noiseZ, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
    fprintf(stdout,"Student-t Noise evidence %lf\n", noiseZ);
165 166 167 168

   } else if (LALInferenceGetProcParamVal(commandLine, "--margphi")) {
    fprintf(stderr, "Using marginalised phase likelihood.\n");
    runState->likelihood=&LALInferenceMarginalisedPhaseLogLikelihood;
169 170 171
   } else if (LALInferenceGetProcParamVal(commandLine, "--margtime")) {
    fprintf(stderr, "Using marginalised time likelihood.\n");
    runState->likelihood=&LALInferenceMarginalisedTimeLogLikelihood;
172 173
   } else if (LALInferenceGetProcParamVal(commandLine, "--margtimephi")) {
     fprintf(stderr, "Using marginalised in time and phase likelihood.\n");
174 175
     runState->likelihood=&LALInferenceMarginalisedTimePhaseLogLikelihood;
     //LALInferenceAddVariable(runState->currentParams, "margtimephi", &margphi, LALINFERENCE_UINT4_t,LALINFERENCE_PARAM_FIXED);
176
   }
177 178 179 180 181
     else if (LALInferenceGetProcParamVal(commandLine, "--roqtime_steps")) {
     fprintf(stderr, "Using ROQ in likelihood.\n");
     runState->likelihood=&LALInferenceUndecomposedFreqDomainLogLikelihood;
    }
     else if (LALInferenceGetProcParamVal(commandLine, "--fastSineGaussianLikelihood")){
182
      fprintf(stderr, "WARNING: Using Fast SineGaussian likelihood and WF for LIB.\n");
183
      runState->likelihood=&LALInferenceFastSineGaussianLogLikelihood;
184
   } else {
Ben Farr's avatar
Ben Farr committed
185
      runState->likelihood=&LALInferenceUndecomposedFreqDomainLogLikelihood;
186 187
   }

188
   /* Try to determine a model-less likelihood, if such a thing makes sense */
John Douglas Veitch's avatar
John Douglas Veitch committed
189
   if (runState->likelihood==&LALInferenceUndecomposedFreqDomainLogLikelihood || runState->likelihood==&LALInferenceMarginalisedPhaseLogLikelihood ){
190 191

                nullLikelihood = LALInferenceNullLogLikelihood(runState->data);
192

193 194
    }
    else if (runState->likelihood==&LALInferenceFreqDomainStudentTLogLikelihood ||
Ben Farr's avatar
Ben Farr committed
195
       runState->likelihood==&LALInferenceMarginalisedTimeLogLikelihood ||
John Douglas Veitch's avatar
John Douglas Veitch committed
196
       runState->likelihood==&LALInferenceMarginalisedTimePhaseLogLikelihood) {
Ben Farr's avatar
Ben Farr committed
197

198
       void *oldtemplate = runState->threads[0]->model->templt;
199 200 201 202 203
       if (runState->threads[0]->model->domain == LAL_SIM_DOMAIN_FREQUENCY) {
         runState->threads[0]->model->templt = &LALInferenceTemplateNullFreqdomain;
       } else {
         runState->threads[0]->model->templt = &LALInferenceTemplateNullTimedomain;
       }
204 205
       nullLikelihood = runState->likelihood(thread->currentParams, runState->data, thread->model);
       runState->threads[0]->model->templt = oldtemplate;
206

Ben Farr's avatar
Ben Farr committed
207
   }
208 209

    //null log likelihood logic doesn't work with noise parameters
Ben Farr's avatar
Ben Farr committed
210
    if (LALInferenceGetProcParamVal(runState->commandLine,"--psdFit") ||
Ben Farr's avatar
Ben Farr committed
211 212 213
       LALInferenceGetProcParamVal(runState->commandLine,"--psd-fit") ||
       LALInferenceGetProcParamVal(runState->commandLine,"--glitchFit") ||
       LALInferenceGetProcParamVal(runState->commandLine,"--glitch-fit")) {
214 215 216 217 218 219
           nullLikelihood = 0.0;
           ifo = runState->data;
           while (ifo != NULL) {
               ifo->nullloglikelihood = 0.0;
               ifo = ifo->next;
           }
Ben Farr's avatar
Ben Farr committed
220
    }
221

222 223
   INT4 t;
   for(t=0; t < runState->nthreads; t++)
224
       runState->threads[t]->nullLikelihood = nullLikelihood;
225

226
   LALInferenceAddVariable(runState->proposalArgs, "nullLikelihood", &nullLikelihood,
227
                           LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
228

229 230 231
    return;
}

232

233 234
const char *non_intrinsic_params[] = {"rightascension", "declination", "polarisation", "time",
                                "deltaLogL", "logL", "deltaloglH1", "deltaloglL1", "deltaloglV1",
235
                                "logw", "logPrior","hrss","loghrss", NULL};
236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260

LALInferenceVariables LALInferenceGetInstrinsicParams(LALInferenceVariables *currentParams)
/***************************************************************/
/* Return a variables structure containing only intrinsic      */
/* parameters.                                                 */
/***************************************************************/
{
    // TODO: add pointer to template function here.
    // (otherwise same parameters but different template will lead to no re-computation!!)
    LALInferenceVariables intrinsicParams;
    const char **non_intrinsic_param = non_intrinsic_params;

    intrinsicParams.head      = NULL;
    intrinsicParams.dimension = 0;
    LALInferenceCopyVariables(currentParams, &intrinsicParams);

    while (*non_intrinsic_param) {
        if (LALInferenceCheckVariable(&intrinsicParams, *non_intrinsic_param))
            LALInferenceRemoveVariable(&intrinsicParams, *non_intrinsic_param);
        non_intrinsic_param++;
    }

    return intrinsicParams;
}

261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276
/* Check to see if item is in the NULL-terminated array.
 If so, return 1. Otherwise, add it to the array and return 0
 */
static int checkItemAndAdd(void *item, void **array);
static int checkItemAndAdd(void *item, void **array)
{
  UINT4 i=0;
  if(!array || !item) return 0;
  while(array[i])
  {
    if(array[i++]==item) return 1;
  }
  array[i]=item;
  return 0;
}

277
/* ============ Likelihood computations: ========== */
278

279 280 281 282
/**
 * For testing purposes (for instance sampling the prior), likelihood that returns 0.0 = log(1) every
 * time.  Activated with the --zeroLogLike command flag.
 */
283
REAL8 LALInferenceZeroLogLikelihood(LALInferenceVariables *currentParams,
284 285
                                    LALInferenceIFOData UNUSED *data,
                                    LALInferenceModel UNUSED *model) {
286 287 288 289 290 291

    INT4 SKY_FRAME=0;
    REAL8 ra,dec,GPSdouble;

    if(LALInferenceCheckVariable(currentParams,"SKY_FRAME"))
      SKY_FRAME=*(INT4 *)LALInferenceGetVariable(currentParams,"SKY_FRAME");
292

293 294 295 296 297 298 299 300 301 302 303
    if(SKY_FRAME==1)
    {
      REAL8 t0=LALInferenceGetREAL8Variable(currentParams,"t0");
      REAL8 alph=acos(LALInferenceGetREAL8Variable(currentParams,"cosalpha"));
      REAL8 theta=LALInferenceGetREAL8Variable(currentParams,"azimuth");
      LALInferenceDetFrameToEquatorial(data->detector,data->next->detector,
                                       t0,alph,theta,&GPSdouble,&ra,&dec);
      LALInferenceAddVariable(currentParams,"rightascension",&ra,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
      LALInferenceAddVariable(currentParams,"declination",&dec,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
      LALInferenceAddVariable(currentParams,"time",&GPSdouble,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
    }
304 305 306 307

  REAL8 zero = 0.0;
  LALInferenceAddVariable(currentParams,"optimal_snr",&zero,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
  LALInferenceAddVariable(currentParams,"matched_filter_snr",&zero,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
308 309 310
  return 0.0;
}

311

312

313

314
REAL8 LALInferenceUndecomposedFreqDomainLogLikelihood(LALInferenceVariables *currentParams,
315 316 317 318 319 320 321 322 323 324 325
                                                      LALInferenceIFOData *data,
                                                      LALInferenceModel *model)
{
  return LALInferenceFusedFreqDomainLogLikelihood(currentParams,
                                                 data,
                                                 model,
                                                  GAUSSIAN);
}


static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *currentParams,
326
                                                        LALInferenceIFOData *data,
327 328
                                                        LALInferenceModel *model,
                                                        LALInferenceLikelihoodFlags marginalisationflags)
329 330 331 332 333 334 335 336 337 338 339 340
/***************************************************************/
/* (log-) likelihood function.                                 */
/* Returns the non-normalised logarithmic likelihood.          */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/* Required (`currentParams') parameters are:                  */
/*   - "rightascension"  (REAL8, radian, 0 <= RA <= 2pi)       */
/*   - "declination"     (REAL8, radian, -pi/2 <= dec <=pi/2)  */
/*   - "polarisation"    (REAL8, radian, 0 <= psi <= ?)        */
/*   - "time"            (REAL8, GPS sec.)                     */
/***************************************************************/
{
  double Fplus, Fcross;
341 342
  //double diffRe, diffIm;
  //double dataReal, dataImag;
343
  double glitchReal=0.0, glitchImag=0.0;
344 345 346
  //REAL8 plainTemplateReal, plainTemplateImag;
  //REAL8 templateReal=0.0, templateImag=0.0;
  int i, j, lower, upper, ifo;
347
  LALInferenceIFOData *dataPtr;
348
  double ra=0.0, dec=0.0, psi=0.0, gmst=0.0;
349
  double GPSdouble=0.0, t0=0.0;
350 351 352
  LIGOTimeGPS GPSlal;
  double chisquared;
  double timedelay;  /* time delay b/w iterferometer & geocenter w.r.t. sky location */
353
  double timeshift=0;  /* time shift (not necessarily same as above)                   */
354
  double deltaT, TwoDeltaToverN, deltaF, twopit=0.0, re, im, dre, dim, newRe, newIm;
355
  double timeTmp;
356
  double mc;
357 358
  /* Burst templates are generated at hrss=1, thus need to rescale amplitude */
  double amp_prefactor=1.0;
359 360 361 362

  COMPLEX16FrequencySeries *calFactor = NULL;
  COMPLEX16 calF = 0.0;

363
  REAL8Vector *logfreqs = NULL;
364 365 366
  REAL8Vector *amps = NULL;
  REAL8Vector *phases = NULL;

367
  UINT4 spcal_active = 0;
368 369 370 371 372
  REAL8 calamp=0.0;
  REAL8 calpha=0.0;
  REAL8 cos_calpha=cos(calpha);
  REAL8 sin_calpha=sin(calpha);
  UINT4 constantcal_active=0;
373
  INT4 errnum=0;
374

375 376
  /* ROQ likelihood stuff */
  REAL8 d_inner_h=0.0;
377
  double dist_min, dist_max;
378
  int cosmology=0;
379 380 381 382 383
  UINT4 margdist = 0;
  if(LALInferenceCheckVariable(model->params, "MARGDIST") && LALInferenceGetVariable(model->params, "MARGDIST"))
  {
      margdist = 1;
      LALInferenceGetMinMaxPrior(model->params, "logdistance", &dist_min, &dist_max);
John Douglas Veitch's avatar
John Douglas Veitch committed
384 385
      dist_max=exp(dist_max);
      dist_min=exp(dist_min);
386 387 388
      if(LALInferenceCheckVariable(model->params,"MARGDIST_COSMOLOGY"))
        cosmology=LALInferenceGetINT4Variable(currentParams,"MARGDIST_COSMOLOGY");

389
  }
390

391 392 393
  if (LALInferenceCheckVariable(currentParams, "spcal_active") && (*(UINT4 *)LALInferenceGetVariable(currentParams, "spcal_active"))) {
    spcal_active = 1;
  }
394 395 396 397 398 399 400
  if (LALInferenceCheckVariable(currentParams, "constantcal_active") && (*(UINT4 *)LALInferenceGetVariable(currentParams, "constantcal_active"))) {
   constantcal_active = 1;
  }
  if (spcal_active && constantcal_active){
    fprintf(stderr,"ERROR: cannot use spline and constant calibration error marginalization together. Exiting...\n");
    exit(1);
  }
401 402 403 404 405
  if (model->roq_flag && constantcal_active){
    fprintf(stderr,"ERROR: cannot use ROQ likelihood and constant calibration error marginalization together. Exiting...\n");
    exit(1);
  }

406 407 408 409 410 411 412 413
  REAL8 degreesOfFreedom=2.0;
  REAL8 chisq=0.0;
  /* margphi params */
  //REAL8 Rre=0.0,Rim=0.0;
  REAL8 D=0.0,S=0.0;
  COMPLEX16 Rcplx=0.0;
  int margphi=0;
  int margtime=0;
414
  REAL8 desired_tc=0.0;
415 416 417 418
  if (marginalisationflags==MARGPHI || marginalisationflags==MARGTIMEPHI)
    margphi=1;
  if (marginalisationflags==MARGTIME || marginalisationflags==MARGTIMEPHI)
    margtime=1;
419

John Douglas Veitch's avatar
John Douglas Veitch committed
420 421 422
  if(model->roq_flag && margtime) XLAL_ERROR_REAL8(XLAL_EINVAL,"ROQ does not support time marginalisation");

  
423 424 425
  LALStatus status;
  memset(&status,0,sizeof(status));

426 427
  if(data==NULL) {XLAL_ERROR_REAL8(XLAL_EINVAL,"ERROR: Encountered NULL data pointer in likelihood\n");}

428 429
  int Nifos=0;
  for(dataPtr=data;dataPtr;dataPtr=dataPtr->next) Nifos++;
430
  void *generatedFreqModels[1+Nifos];
431
  for(i=0;i<=Nifos;i++) generatedFreqModels[i]=NULL;
432

433 434
  //noise model meta parameters
  gsl_matrix *nparams = NULL;//pointer to matrix holding noise parameters
435 436 437 438

  gsl_matrix *psdBandsMin  = NULL;//pointer to matrix holding min frequencies for psd model
  gsl_matrix *psdBandsMax = NULL;//pointer to matrix holding max frequencies for psd model

439
  //different formats for storing glitch model for DWT, FFT, and integration
440
  gsl_matrix *glitchFD=NULL;
441

442
  int Nblock = 1;            //number of frequency blocks per IFO
443
  int psdFlag = 0;           //flag for including psd fitting
444 445
  int glitchFlag = 0;   //flag for including glitch model
  int signalFlag = 1;   //flag for including signal model
446 447

  //check if psd parameters are included in the model
448
  psdFlag = 0;
449 450
  if(LALInferenceCheckVariable(currentParams, "psdScaleFlag"))
    psdFlag = *((INT4 *)LALInferenceGetVariable(currentParams, "psdScaleFlag"));
451 452 453 454 455 456
  if(psdFlag)
  {
    //if so, store current noise parameters in easily accessible matrix
    nparams = *((gsl_matrix **)LALInferenceGetVariable(currentParams, "psdscale"));
    Nblock = (int)nparams->size2;

457 458 459
    psdBandsMin = *((gsl_matrix **)LALInferenceGetVariable(currentParams, "psdBandsMin"));
    psdBandsMax = *((gsl_matrix **)LALInferenceGetVariable(currentParams, "psdBandsMax"));

460 461 462 463
  }
  double alpha[Nblock];
  double lnalpha[Nblock];

464 465 466
  double psdBandsMin_array[Nblock];
  double psdBandsMax_array[Nblock];

467 468 469 470 471 472 473 474 475 476 477
  //check if glitch model is being used
  glitchFlag = 0;
  if(LALInferenceCheckVariable(currentParams,"glitchFitFlag"))
    glitchFlag = *((INT4 *)LALInferenceGetVariable(currentParams, "glitchFitFlag"));
  if(glitchFlag)
    glitchFD = *((gsl_matrix **)LALInferenceGetVariable(currentParams, "morlet_FD"));

  //check if signal model is being used
  signalFlag=1;
  if(LALInferenceCheckVariable(currentParams, "signalModelFlag"))
    signalFlag = *((INT4 *)LALInferenceGetVariable(currentParams, "signalModelFlag"));
478

479 480 481 482 483 484 485 486 487 488 489 490 491 492 493
  int freq_length=0,time_length=0;
  COMPLEX16Vector * dh_S_tilde=NULL;
  COMPLEX16Vector * dh_S_phase_tilde = NULL;
  REAL8Vector *dh_S=NULL;
  REAL8Vector *dh_S_phase = NULL;
  /* Setup times to integrate over */
  freq_length = data->freqData->data->length;
  time_length = 2*(freq_length-1);

  /* Desired tc == 2 seconds before buffer end.  Only used during
     margtime{phi} to try to place the waveform in a reasonable
     place before time-shifting */
  deltaT = data->timeData->deltaT;
  REAL8 epoch = XLALGPSGetREAL8(&(data->freqData->epoch));
  desired_tc = epoch + (time_length-1)*deltaT - 2.0;
494

495 496
  if(signalFlag)
  {
497 498 499 500
    if(LALInferenceCheckVariable(currentParams,"logmc")){
      mc=exp(*(REAL8 *)LALInferenceGetVariable(currentParams,"logmc"));
      LALInferenceAddVariable(currentParams,"chirpmass",&mc,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
    }
501 502 503 504 505 506
    if(LALInferenceCheckVariable(currentParams, "loghrss")){
      amp_prefactor = exp(*(REAL8*)LALInferenceGetVariable(currentParams,"loghrss"));
    }
    else if (LALInferenceCheckVariable(currentParams, "hrss")){
      amp_prefactor = (*(REAL8*)LALInferenceGetVariable(currentParams,"hrss"));
    }
507

508 509 510 511 512 513 514 515
    INT4 SKY_FRAME=0;
    if(LALInferenceCheckVariable(currentParams,"SKY_FRAME"))
      SKY_FRAME=*(INT4 *)LALInferenceGetVariable(currentParams,"SKY_FRAME");
    if(SKY_FRAME==0){
      /* determine source's sky location & orientation parameters: */
      ra        = *(REAL8*) LALInferenceGetVariable(currentParams, "rightascension"); /* radian      */
      dec       = *(REAL8*) LALInferenceGetVariable(currentParams, "declination");    /* radian      */
    }
516
    else
517
    {
518 519 520 521
	    if(Nifos<2){
		    fprintf(stderr,"ERROR: Cannot use --detector-frame with less than 2 detectors!\n");
		    exit(1);
	    }
522 523 524 525 526 527 528 529
      if(!margtime)
      {
         t0=LALInferenceGetREAL8Variable(currentParams,"t0");
      }
      else /* Use the desired end time to compute the mapping to ra,dec */
      {
              t0=desired_tc;
      }
530 531 532 533 534 535
      REAL8 alph=acos(LALInferenceGetREAL8Variable(currentParams,"cosalpha"));
      REAL8 theta=LALInferenceGetREAL8Variable(currentParams,"azimuth");
      LALInferenceDetFrameToEquatorial(data->detector,data->next->detector,
                                       t0,alph,theta,&GPSdouble,&ra,&dec);
      LALInferenceAddVariable(currentParams,"rightascension",&ra,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
      LALInferenceAddVariable(currentParams,"declination",&dec,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
536
      if(!margtime) LALInferenceAddVariable(currentParams,"time",&GPSdouble,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
537 538
    }
    psi       = *(REAL8*) LALInferenceGetVariable(currentParams, "polarisation");   /* radian      */
539 540 541 542
    if(!margtime)
	      GPSdouble = *(REAL8*) LALInferenceGetVariable(currentParams, "time");           /* GPS seconds */
    else
	      GPSdouble = XLALGPSGetREAL8(&(data->freqData->epoch));
543

544

545 546 547 548 549 550 551
    // Add phase parameter set to 0 for calculation
    if(margphi ){
      REAL8 phi0=0.0;
      if(LALInferenceCheckVariable(currentParams,"phase")) LALInferenceRemoveVariable(currentParams,"phase");
      LALInferenceAddVariable(currentParams, "phase",&phi0,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
    }
  }
552

553

554 555
  if(margtime)
  {
556
    GPSdouble = desired_tc;
557 558
    dh_S_tilde = XLALCreateCOMPLEX16Vector(freq_length);
    dh_S = XLALCreateREAL8Vector(time_length);
559

560 561
    if (dh_S_tilde ==NULL || dh_S == NULL)
      XLAL_ERROR_REAL8(XLAL_ENOMEM, "Out of memory in LALInferenceMarginalisedTimeLogLikelihood.");
562

563
    for (i = 0; i < freq_length; i++) {
564 565
      dh_S_tilde->data[i] = 0.0;
    }
566 567 568 569 570 571 572 573 574 575 576 577 578

    if (margphi) {
      dh_S_phase_tilde = XLALCreateCOMPLEX16Vector(freq_length);
      dh_S_phase = XLALCreateREAL8Vector(time_length);

      if (dh_S_phase_tilde == NULL || dh_S_phase == NULL) {
	XLAL_ERROR_REAL8(XLAL_ENOMEM, "Out of memory in time-phase marginalised likelihood.");
      }

      for (i = 0; i < freq_length; i++) {
	dh_S_phase_tilde->data[i] = 0.0;
      }
    }
579
  }
580

581 582 583
  /* figure out GMST: */
  XLALGPSSetREAL8(&GPSlal, GPSdouble);
  gmst=XLALGreenwichMeanSiderealTime(&GPSlal);
584

585
  chisquared = 0.0;
586 587 588 589 590
  REAL8 loglikelihood = 0.0;

  /* Reset SNR */
  model->SNR = 0.0;

591
  /* loop over data (different interferometers): */
592
  for(dataPtr=data,ifo=0; dataPtr; dataPtr=dataPtr->next,ifo++) {
593 594 595 596 597
    /* The parameters the Likelihood function can handle by itself   */
    /* (and which shouldn't affect the template function) are        */
    /* sky location (ra, dec), polarisation and signal arrival time. */
    /* Note that the template function shifts the waveform to so that*/
	/* t_c corresponds to the "time" parameter in                    */
598
	/* model->params (set, e.g., from the trigger value).     */
599

600
    /* Reset log-likelihood */
601 602
    model->ifo_loglikelihoods[ifo] = 0.0;
    model->ifo_SNRs[ifo] = 0.0;
603 604
    COMPLEX16 this_ifo_d_inner_h = 0.0;
    REAL8 this_ifo_s = 0.0;
605 606 607 608
    // Check if student-t likelihood is being used
    if(marginalisationflags==STUDENTT)
    {
      /* extract the element from the "df" vector that carries the current Ifo's name: */
609
      CHAR df_variable_name[320];
610 611 612 613 614 615 616 617 618 619 620 621 622
      snprintf(df_variable_name,sizeof(df_variable_name),"df_%s",dataPtr->name);
      if(LALInferenceCheckVariable(currentParams,df_variable_name)){
        printf("Found variable %s\n",df_variable_name);
        degreesOfFreedom = *(REAL8*) LALInferenceGetVariable(currentParams,df_variable_name);
      }
      else {
        degreesOfFreedom = dataPtr->STDOF;
      }
      if (!(degreesOfFreedom>0)) {
        XLALPrintError(" ERROR in StudentTLogLikelihood(): degrees-of-freedom parameter must be positive.\n");
        XLAL_ERROR_REAL8(XLAL_EDOM);
      }
    }
623

624
    if(signalFlag){
625

626 627 628 629 630 631 632 633 634 635 636 637
      /* Check to see if this buffer has already been filled with the signal.
       Different dataPtrs can share the same signal buffer to avoid repeated
       calls to template */
      if(!checkItemAndAdd((void *)(model->freqhPlus), generatedFreqModels))
      {
        /* Compare parameter values with parameter values corresponding  */
        /* to currently stored template; ignore "time" variable:         */
        if (LALInferenceCheckVariable(model->params, "time")) {
          timeTmp = *(REAL8 *) LALInferenceGetVariable(model->params, "time");
          LALInferenceRemoveVariable(model->params, "time");
        }
        else timeTmp = GPSdouble;
638

639 640 641 642
        LALInferenceCopyVariables(currentParams, model->params);
        // Remove time variable so it can be over-written (if it was pinned)
        if(LALInferenceCheckVariable(model->params,"time")) LALInferenceRemoveVariable(model->params,"time");
        LALInferenceAddVariable(model->params, "time", &timeTmp, LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_LINEAR);
643

644 645 646 647 648 649
        XLAL_TRY(model->templt(model),errnum);
        errnum&=~XLAL_EFUNC;
        if(errnum!=XLAL_SUCCESS)
        {
          switch(errnum)
          {
650
            case XLAL_EUSR0: /* Template generation failed in a known way, set -Inf likelihood */
651 652 653 654 655
		      /* Free up allocated vectors */
              if(dh_S_tilde) XLALDestroyCOMPLEX16Vector(dh_S_tilde);
              if(dh_S) XLALDestroyREAL8Vector(dh_S);
              if(dh_S_phase_tilde) XLALDestroyCOMPLEX16Vector(dh_S_phase_tilde);
              if(dh_S_phase) XLALDestroyREAL8Vector(dh_S_phase);
656 657 658 659 660 661 662
              if(model->roq_flag)
              {
                if ( model->roq->hptildeLinear ) XLALDestroyCOMPLEX16FrequencySeries(model->roq->hptildeLinear);
                if ( model->roq->hctildeLinear ) XLALDestroyCOMPLEX16FrequencySeries(model->roq->hctildeLinear);
                if ( model->roq->hptildeQuadratic ) XLALDestroyCOMPLEX16FrequencySeries(model->roq->hptildeQuadratic);
                if ( model->roq->hctildeQuadratic ) XLALDestroyCOMPLEX16FrequencySeries(model->roq->hctildeQuadratic);
              }
663
              return (-INFINITY);
664 665 666 667 668 669 670
              break;
            default: /* Panic! */
              fprintf(stderr,"Unhandled error in template generation - exiting!\n");
              fprintf(stderr,"XLALError: %d, %s\n",errnum,XLALErrorString(errnum));
              exit(1);
              break;
          }
671

672
        }
673

674 675 676 677 678
        if (model->domain == LAL_SIM_DOMAIN_TIME) {
          /* TD --> FD. */
          LALInferenceExecuteFT(model);
        }
      }
679

680
        /* Template is now in model->timeFreqhPlus and hCross */
681

682 683 684
        /* Calibration stuff if necessary */
        /*spline*/
        if (spcal_active) {
685
          logfreqs = NULL;
686
          amps = NULL;
687 688 689 690
          phases = NULL;
	  /* get_calib_spline creates and fills the logfreqs, amps, phases arrays */
	  get_calib_spline(currentParams, dataPtr->name, &logfreqs, &amps, &phases);
	  if (model->roq_flag) {
691 692 693

             LALInferenceSplineCalibrationFactorROQ(logfreqs, amps, phases,
						model->roq->frequencyNodesLinear,
694
						&(model->roq->calFactorLinear),
695
						model->roq->frequencyNodesQuadratic,
696
						&(model->roq->calFactorQuadratic));
697
	  }
698

699 700 701
	  else{
	    if (calFactor == NULL) {
	      calFactor = XLALCreateCOMPLEX16FrequencySeries("calibration factors",
702 703 704 705
                       &(dataPtr->freqData->epoch),
                       0, dataPtr->freqData->deltaF,
                       &lalDimensionlessUnit,
                       dataPtr->freqData->data->length);
706
	    }
707
          LALInferenceSplineCalibrationFactor(logfreqs, amps, phases, calFactor);
708
	}
709 710 711
	if(logfreqs) XLALDestroyREAL8Vector(logfreqs);
	if(amps) XLALDestroyREAL8Vector(amps);
	if(phases) XLALDestroyREAL8Vector(phases);
712

713 714 715
        }
        /*constant*/
        if (constantcal_active){
716
          char CA_A[320];
717 718 719 720
          sprintf(CA_A,"%s_%s","calamp",dataPtr->name);
          if (LALInferenceCheckVariable(currentParams, CA_A))
            calamp=(*(REAL8*) LALInferenceGetVariable(currentParams, CA_A));
          else
721
            calamp=0.0;
722
          char CP_A[320];
723 724 725 726
          sprintf(CP_A,"%s_%s","calpha",dataPtr->name);
          if (LALInferenceCheckVariable(currentParams, CP_A))
            calpha=(*(REAL8*) LALInferenceGetVariable(currentParams, CP_A));
          else
727
            calpha=0.0;
728 729 730
          cos_calpha=cos(calpha);
          sin_calpha=-sin(calpha);
        }
731 732
        /* determine beam pattern response (F_plus and F_cross) for given Ifo: */
        XLALComputeDetAMResponse(&Fplus, &Fcross, (const REAL4(*)[3])dataPtr->detector->response, ra, dec, psi, gmst);
733

734 735 736 737
        /* signal arrival time (relative to geocenter); */
        timedelay = XLALTimeDelayFromEarthCenter(dataPtr->detector->location, ra, dec, &GPSlal);
        /* (negative timedelay means signal arrives earlier at Ifo than at geocenter, etc.) */
        /* amount by which to time-shift template (not necessarily same as above "timedelay"): */
738 739 740 741 742 743 744
        if (margtime)
          /* If we are marginalising over time, we want the
	      freq-domain signal to have tC = epoch, so we shift it
	      from the model's "time" parameter to epoch */
          timeshift =  (epoch - (*(REAL8 *) LALInferenceGetVariable(model->params, "time"))) + timedelay;
        else
          timeshift =  (GPSdouble - (*(REAL8*) LALInferenceGetVariable(model->params, "time"))) + timedelay;
745
        twopit    = LAL_TWOPI * timeshift;
746

747 748 749 750
        /* For burst, add the right hrss in the amplitude. */
        Fplus*=amp_prefactor;
        Fcross*=amp_prefactor;

751 752 753
        dataPtr->fPlus = Fplus;
        dataPtr->fCross = Fcross;
        dataPtr->timeshift = timeshift;
754
    }//end signalFlag condition
755 756 757 758

    /* determine frequency range & loop over frequency bins: */
    deltaT = dataPtr->timeData->deltaT;
    deltaF = 1.0 / (((double)dataPtr->timeData->data->length) * deltaT);
759 760
    lower = (UINT4)ceil(dataPtr->fLow / deltaF);
    upper = (UINT4)floor(dataPtr->fHigh / deltaF);
761
    TwoDeltaToverN = 2.0 * deltaT / ((double) dataPtr->timeData->data->length);
762 763 764 765 766 767 768 769

    /* Employ a trick here for avoiding cos(...) and sin(...) in time
       shifting.  We need to multiply each template frequency bin by
       exp(-J*twopit*deltaF*i) = exp(-J*twopit*deltaF*(i-1)) +
       exp(-J*twopit*deltaF*(i-1))*(exp(-J*twopit*deltaF) - 1) .  This
       recurrance relation has the advantage that the error growth is
       O(sqrt(N)) for N repetitions. */

Vivien Raymond's avatar
Vivien Raymond committed
770
    /* See, for example,
771 772

       Press, Teukolsky, Vetteling & Flannery, 2007.  Numerical
Vivien Raymond's avatar
Vivien Raymond committed
773
       Recipes, Third Edition, Chapter 5.4.
774 775 776

       Singleton, 1967. On computing the fast Fourier
       transform. Comm. ACM, vol. 10, 647–654. */
Vivien Raymond's avatar
Vivien Raymond committed
777

778 779 780
    /* Incremental values, using cos(theta) - 1 = -2*sin(theta/2)^2 */
    dim = -sin(twopit*deltaF);
    dre = -2.0*sin(0.5*twopit*deltaF)*sin(0.5*twopit*deltaF);
781

782
    //Set up noise PSD meta parameters
783
    for(i=0; i<Nblock; i++)
784
    {
785
      if(psdFlag)
786 787 788
      {
        alpha[i]   = gsl_matrix_get(nparams,ifo,i);
        lnalpha[i] = log(alpha[i]);
789 790 791

        psdBandsMin_array[i] = gsl_matrix_get(psdBandsMin,ifo,i);
        psdBandsMax_array[i] = gsl_matrix_get(psdBandsMax,ifo,i);
792
      }
793 794 795 796 797
      else
      {
        alpha[i]=1.0;
        lnalpha[i]=0.0;
      }
798 799
    }

800 801
    if (model->roq_flag) {

802
	double complex weight_iii;
803 804

	if (spcal_active){
805

806
	    for(unsigned int iii=0; iii < model->roq->frequencyNodesLinear->length; iii++){
807

808 809 810 811 812
			complex double template_EI = model->roq->calFactorLinear->data[iii] * (dataPtr->fPlus*model->roq->hptildeLinear->data->data[iii] + dataPtr->fCross*model->roq->hctildeLinear->data->data[iii] );

			weight_iii = gsl_spline_eval (dataPtr->roq->weights_linear[iii].spline_real_weight_linear, timeshift, dataPtr->roq->weights_linear[iii].acc_real_weight_linear) + I*gsl_spline_eval (dataPtr->roq->weights_linear[iii].spline_imag_weight_linear, timeshift, dataPtr->roq->weights_linear[iii].acc_imag_weight_linear);

			this_ifo_d_inner_h += ( weight_iii * ( conj( template_EI ) ) );
813
		}
814

815 816 817 818 819
		for(unsigned int jjj=0; jjj < model->roq->frequencyNodesQuadratic->length; jjj++){

			this_ifo_s += dataPtr->roq->weightsQuadratic[jjj] * creal( conj( model->roq->calFactorQuadratic->data[jjj] * (model->roq->hptildeQuadratic->data->data[jjj]*dataPtr->fPlus + model->roq->hctildeQuadratic->data->data[jjj]*dataPtr->fCross) ) * ( model->roq->calFactorQuadratic->data[jjj] * (model->roq->hptildeQuadratic->data->data[jjj]*dataPtr->fPlus + model->roq->hctildeQuadratic->data->data[jjj]*dataPtr->fCross) ) );
		}
	}
820

821
	else{
822

823 824
		for(unsigned int iii=0; iii < model->roq->frequencyNodesLinear->length; iii++){

825 826
			complex double template_EI = dataPtr->fPlus*model->roq->hptildeLinear->data->data[iii] + dataPtr->fCross*model->roq->hctildeLinear->data->data[iii];

827
			weight_iii = gsl_spline_eval (dataPtr->roq->weights_linear[iii].spline_real_weight_linear, timeshift, dataPtr->roq->weights_linear[iii].acc_real_weight_linear) + I*gsl_spline_eval (dataPtr->roq->weights_linear[iii].spline_imag_weight_linear, timeshift, dataPtr->roq->weights_linear[iii].acc_imag_weight_linear);
828 829 830

			this_ifo_d_inner_h += weight_iii*conj(template_EI) ;

831
		}
832

833
		for(unsigned int jjj=0; jjj < model->roq->frequencyNodesQuadratic->length; jjj++){
834
			complex double template_EI = model->roq->hptildeQuadratic->data->data[jjj]*Fplus + model->roq->hctildeQuadratic->data->data[jjj]*Fcross;
835

836 837
			this_ifo_s += dataPtr->roq->weightsQuadratic[jjj] * creal( conj(template_EI) * (template_EI) );
					}
838 839
	}

840 841 842 843 844 845 846 847
    d_inner_h += creal(this_ifo_d_inner_h);
    // D gets the factor of 2 inside nullloglikelihood
    // R gets a factor of 2 later, before entering the Bessel function
    S += 0.5*this_ifo_s;
    D += -dataPtr->nullloglikelihood;
    Rcplx += 0.5*this_ifo_d_inner_h;

    model->ifo_loglikelihoods[ifo] = creal(this_ifo_d_inner_h) - (0.5*this_ifo_s) + dataPtr->nullloglikelihood;
John Douglas Veitch's avatar
John Douglas Veitch committed
848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865
    
    switch(marginalisationflags)
    {
    case GAUSSIAN:
    case STUDENTT:
      loglikelihood += model->ifo_loglikelihoods[ifo];
      break;
    case MARGTIME:
    case MARGPHI:
    case MARGTIMEPHI:
      /* These are non-separable likelihoods, so single IFO log(L)
	 doesn't make sense. */
      model->ifo_loglikelihoods[ifo] = 0.0;
      break;
    default:
      break;
    }
    
866
	char varname[VARNAME_MAX];
867 868 869 870 871 872 873
    if((VARNAME_MAX <= snprintf(varname,VARNAME_MAX,"%s_optimal_snr",dataPtr->name)))
    {
        fprintf(stderr,"variable name too long\n"); exit(1);
    }
    REAL8 this_ifo_snr = sqrt(this_ifo_s);
    model->ifo_SNRs[ifo] = this_ifo_snr;
    LALInferenceAddREAL8Variable(currentParams,varname,this_ifo_snr,LALINFERENCE_PARAM_OUTPUT);
874

875 876 877 878 879 880
    if((VARNAME_MAX <= snprintf(varname,VARNAME_MAX,"%s_cplx_snr_amp",dataPtr->name)))
    {
        fprintf(stderr,"variable name too long\n"); exit(1);
    }
    REAL8 cplx_snr_amp = cabs(this_ifo_d_inner_h)/this_ifo_snr;
    LALInferenceAddREAL8Variable(currentParams,varname,cplx_snr_amp,LALINFERENCE_PARAM_OUTPUT);
881

882 883 884 885 886 887
    if((VARNAME_MAX <= snprintf(varname,VARNAME_MAX,"%s_cplx_snr_arg",dataPtr->name)))
    {
        fprintf(stderr,"variable name too long\n"); exit(1);
    }
    REAL8 cplx_snr_phase = carg(this_ifo_d_inner_h);
    LALInferenceAddREAL8Variable(currentParams,varname,cplx_snr_phase,LALINFERENCE_PARAM_OUTPUT);
888

889
    }
890
    else{
891
    REAL8 *psd=&(dataPtr->oneSidedNoisePowerSpectrum->data->data[lower]);
892 893 894 895 896 897
    COMPLEX16 *dtilde=&(dataPtr->freqData->data->data[lower]);
    COMPLEX16 *hptilde=&(model->freqhPlus->data->data[lower]);
    COMPLEX16 *hctilde=&(model->freqhCross->data->data[lower]);
    COMPLEX16 diff=0.0;
    COMPLEX16 template=0.0;
    REAL8 templatesq=0.0;
898
    REAL8 this_ifo_S=0.0;
899
    COMPLEX16 this_ifo_Rcplx=0.0;
900

901 902 903 904 905 906 907
    for (i=lower,chisq=0.0,re = cos(twopit*deltaF*i),im = -sin(twopit*deltaF*i);
         i<=upper;
         i++, psd++, hptilde++, hctilde++, dtilde++,
         newRe = re + re*dre - im*dim,
         newIm = im + re*dim + im*dre,
         re = newRe, im = newIm)
    {
908

909 910 911 912
      COMPLEX16 d=*dtilde;
      /* Normalise PSD to our funny standard (see twoDeltaTOverN
	 below). */
      REAL8 sigmasq=(*psd)*deltaT*deltaT;
913

914
      if (constantcal_active) {
915 916 917 918 919
        REAL8 dre_tmp= creal(d)*cos_calpha - cimag(d)*sin_calpha;
        REAL8 dim_tmp = creal(d)*sin_calpha + cimag(d)*cos_calpha;
        dre_tmp/=(1.0+calamp);
        dim_tmp/=(1.0+calamp);

920 921
        d=crect(dre_tmp,dim_tmp);
        sigmasq/=((1.0+calamp)*(1.0+calamp));
922 923
      }

924
      REAL8 singleFreqBinTerm;
925 926


927 928 929 930 931 932 933 934 935 936 937 938
      /* Add noise PSD parameters to the model */
      if(psdFlag)
      {
        for(j=0; j<Nblock; j++)
        {
          if (i >= psdBandsMin_array[j] && i <= psdBandsMax_array[j])
          {
            sigmasq  *= alpha[j];
            loglikelihood -= lnalpha[j];
          }
        }
      }
939

940
      //subtract GW model from residual
Ben Farr's avatar
Ben Farr committed
941 942
      diff = d;

943
      if(signalFlag){
944
      /* derive template (involving location/orientation parameters) from given plus/cross waveforms: */
945
      COMPLEX16 plainTemplate = Fplus*(*hptilde)+Fcross*(*hctilde);
946

947
      /* Do time shifting */
948
      template = plainTemplate * (re + I*im);
949

950
      if (spcal_active) {
951 952
          calF = calFactor->data->data[i];
          template = template*calF;
953
      }
954

Ben Farr's avatar
Ben Farr committed
955
      diff -= template;
956

957 958 959 960 961 962 963 964
      }//end signal subtraction

      //subtract glitch model from residual
      if(glitchFlag)
      {
        /* fourier amplitudes of glitches */
        glitchReal = gsl_matrix_get(glitchFD,ifo,2*i);
        glitchImag = gsl_matrix_get(glitchFD,ifo,2*i+1);
965
        COMPLEX16 glitch = glitchReal + I*glitchImag;
966
        diff -=glitch*deltaT;
967 968 969

      }//end glitch subtraction

970 971 972
      templatesq=creal(template)*creal(template) + cimag(template)*cimag(template);
      REAL8 datasq = creal(d)*creal(d)+cimag(d)*cimag(d);
      D+=TwoDeltaToverN*datasq/sigmasq;
973
      this_ifo_S+=TwoDeltaToverN*templatesq/sigmasq;
974
      COMPLEX16 dhstar = TwoDeltaToverN*d*conj(template)/sigmasq;
975
      this_ifo_Rcplx+=dhstar;
976
      Rcplx+=dhstar;
Vivien Raymond's avatar
Vivien Raymond committed
977

978
      switch(marginalisationflags)
979
      {
980
        case GAUSSIAN:
981
        {
982 983 984 985 986 987
          REAL8 diffsq = creal(diff)*creal(diff)+cimag(diff)*cimag(diff);
          chisq = TwoDeltaToverN*diffsq/sigmasq;
          singleFreqBinTerm = chisq;
          chisquared  += singleFreqBinTerm;
          model->ifo_loglikelihoods[ifo] -= singleFreqBinTerm;
          break;
988
        }
989 990 991 992 993 994 995 996 997 998 999 1000 1001
        case STUDENTT:
        {
          REAL8 diffsq = creal(diff)*creal(diff)+cimag(diff)*cimag(diff);
          chisq = TwoDeltaToverN*diffsq/sigmasq;
          singleFreqBinTerm = ((degreesOfFreedom+2.0)/2.0) * log(1.0 + chisq/degreesOfFreedom) ;
          chisquared  += singleFreqBinTerm;
          model->ifo_loglikelihoods[ifo] -= singleFreqBinTerm;
          break;
        }
        case MARGTIME:
        case MARGTIMEPHI:
        {
          loglikelihood+=-TwoDeltaToverN*(templatesq+datasq)/sigmasq;
1002

1003 1004 1005 1006 1007 1008
          /* Note: No Factor of 2 here, since we are using the 2-sided
	     COMPLEX16FFT.  Also, we use d*conj(h) because we are
	     using a complex->real *inverse* FFT to compute the
	     time-series of likelihoods. */
          dh_S_tilde->data[i] += TwoDeltaToverN * d * conj(template) / sigmasq;

1009 1010 1011 1012
          if (margphi) {
            /* This is the other phase quadrature */
            dh_S_phase_tilde->data[i] += TwoDeltaToverN * d * conj(I*template) / sigmasq;
          }
1013

1014 1015 1016
          break;
        }
        case MARGPHI:
1017
        {
1018
          break;
1019
        }
1020 1021
        default:
          break;
1022 1023
      }

1024

1025

1026 1027 1028
    } /* End loop over freq bins */
    switch(marginalisationflags)
    {
1029 1030 1031 1032 1033 1034 1035
    case GAUSSIAN:
    case STUDENTT:
      loglikelihood += model->ifo_loglikelihoods[ifo];
      break;
    case MARGTIME:
    case MARGPHI:
    case MARGTIMEPHI:
1036 1037 1038 1039
      /* These are non-separable likelihoods, so single IFO log(L)
	 doesn't make sense. */
      model->ifo_loglikelihoods[ifo] = 0.0;
      break;
1040 1041
    default:
      break;
1042
    }
1043 1044
    S+=this_ifo_S;
    char varname[VARNAME_MAX];
1045 1046 1047 1048
    if((VARNAME_MAX <= snprintf(varname,VARNAME_MAX,"%s_optimal_snr",dataPtr->name)))
    {
        fprintf(stderr,"variable name too long\n"); exit(1);
    }
1049
    LALInferenceAddREAL8Variable(currentParams,varname,sqrt(2.0*this_ifo_S),LALINFERENCE_PARAM_OUTPUT);
1050

1051 1052 1053 1054
    if((VARNAME_MAX <= snprintf(varname,VARNAME_MAX,"%s_cplx_snr_amp",dataPtr->name)))
    {
        fprintf(stderr,"variable name too long\n"); exit(1);
    }
1055 1056 1057 1058 1059
    REAL8 cplx_snr_amp=0.0;
    REAL8 cplx_snr_phase=carg(this_ifo_Rcplx);
    if(this_ifo_S > 0) cplx_snr_amp=2.0*cabs(this_ifo_Rcplx)/sqrt(2.0*this_ifo_S);

    LALInferenceAddREAL8Variable(currentParams,varname,cplx_snr_amp,LALINFERENCE_PARAM_OUTPUT);
1060

1061 1062 1063 1064
    if((VARNAME_MAX <= snprintf(varname,VARNAME_MAX,"%s_cplx_snr_arg",dataPtr->name)))
    {
        fprintf(stderr,"variable name too long\n"); exit(1);
    }
1065
    LALInferenceAddREAL8Variable(currentParams,varname,cplx_snr_phase,LALINFERENCE_PARAM_OUTPUT);
1066 1067
    if(margdist )
      {
1068
          XLAL_TRY(model->ifo_loglikelihoods[ifo] = LALInferenceMarginalDistanceLogLikelihood(dist_min, dist_max, sqrt(this_ifo_S), 2.0*cabs(this_ifo_Rcplx), cosmology), errnum);
1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084
          errnum&=~XLAL_EFUNC;
          if(errnum!=XLAL_SUCCESS)
          {
            switch(errnum)
            {
              case XLAL_ERANGE: /* The SNR input was outside the interpolation range */
		if(calFactor) {XLALDestroyCOMPLEX16FrequencySeries(calFactor); calFactor=NULL;}
                return (-INFINITY);
                break;
              default: /* Panic! */
                fprintf(stderr,"Unhandled error in marginal distance likelihood - exiting!\n");
                fprintf(stderr,"XLALError: %d, %s\n",errnum,XLALErrorString(errnum));
                exit(1);
                break;
            }
          }
1085
      }
1086
   /* Clean up calibration if necessary */
1087 1088 1089 1090
    if (!(calFactor == NULL)) {
      XLALDestroyCOMPLEX16FrequencySeries(calFactor);
      calFactor = NULL;
    }
1091
  } /* end loop over detectors */
1092

1093 1094 1095 1096 1097 1098
  }
  if (model->roq_flag){



	REAL8 OptimalSNR=sqrt(S);
1099
        REAL8 MatchedFilterSNR = d_inner_h/OptimalSNR;