Maintenance will be performed on git.ligo.org, chat.ligo.org, containers.ligo.org, and docs.ligo.org starting 2 March 2020 at approximately 8am MST. It is expected to take around 10 minutes and will include a short period of downtime towards the end of the maintenance window. Please direct any comments, concerns, or questions to computing-help@igwn.org.

Commit c8d9145d authored by John Douglas Veitch's avatar John Douglas Veitch

Fix some normalisation

parent e6c27a12
......@@ -62,6 +62,10 @@ typedef enum
#define UNUSED
#endif
#ifndef _OPENMP
#define omp ignore
#endif
static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *currentParams,
LALInferenceIFOData *data,
LALInferenceModel *model,
......@@ -374,6 +378,8 @@ static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *cur
{
margdist = 1;
LALInferenceGetMinMaxPrior(model->params, "logdistance", &dist_min, &dist_max);
dist_max=exp(dist_max);
dist_min=exp(dist_min);
}
if (LALInferenceCheckVariable(currentParams, "spcal_active") && (*(UINT4 *)LALInferenceGetVariable(currentParams, "spcal_active"))) {
......@@ -1242,31 +1248,31 @@ static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *cur
double LALInferenceMarginalDistanceLogLikelihood(double dist_min, double dist_max, double OptimalSNR, double d_inner_h)
{
static const size_t default_log_radial_integrator_size = 400;
static const size_t default_log_radial_integrator_size = 400;
double loglikelihood=0;
static log_radial_integrator *integrator;
#pragma omp threadprivate(integrator)
double pmax = 100000; /* CHECKME: Max SNR allowed ? */
static log_radial_integrator *integrator;
if (integrator == NULL)
{
double pmax = 100000; /* CHECKME: Max SNR allowed ? */
#pragma omp critical
{
if (integrator == NULL)
{
printf("Initialising distance integration lookup table\n");
int cosmology = 0; /* 0 = euclidean, nonzero co-moving */
/* Initialise the integrator for the first time */
integrator = log_radial_integrator_init(
dist_min,
dist_max,
2, /* Power of distance in prior */
cosmology,
pmax,
default_log_radial_integrator_size * 5 /* CHECKME: fudge factor of 5 compared to bayestar */
);
if (!integrator) XLAL_ERROR(XLAL_EFUNC, "Unable to initialise distance marginalisation integrator");
}
int cosmology = 0; /* 0 = euclidean, nonzero co-moving */
/* Initialise the integrator for the first time */
integrator = log_radial_integrator_init(
dist_min,
dist_max,
2, /* Power of distance in prior */
cosmology,
pmax,
default_log_radial_integrator_size * 5 /* CHECKME: fudge factor of 5 compared to bayestar */
);
}
}
#pragma omp barrier
if (!integrator) XLAL_ERROR(XLAL_EFUNC, "Unable to initialise distance marginalisation integrator");
//double marg_l = dist_integral(OptimalSNR*OptimalSNR, 2.0*d_inner_h, exp(dist_min), exp(dist_max));
if (isnan(OptimalSNR) || isnan(d_inner_h) || pmax<OptimalSNR / sqrt(2))
{
loglikelihood=-INFINITY;
......@@ -1275,7 +1281,7 @@ double LALInferenceMarginalDistanceLogLikelihood(double dist_min, double dist_ma
else
{
double marg_l = log_radial_integrator_eval(integrator, OptimalSNR / sqrt(2), 2.0 * d_inner_h, log(OptimalSNR / sqrt(2)), log(2.0* d_inner_h));
loglikelihood = marg_l;
loglikelihood = marg_l - (log(dist_max-dist_min)); /* Normalise over prior */
}
return (loglikelihood);
}
......
......@@ -216,14 +216,14 @@ double log_radial_integral(double r1, double r2, double p, double b, int k, int
/* Set up integrand data structure. */
const gsl_function func = {radial_integrand, &params};
gsl_error_handler_t *old_handler = gsl_set_error_handler_off();
//gsl_error_handler_t *old_handler = gsl_set_error_handler_off();
/* Perform adaptive Gaussian quadrature. */
ret = gsl_integration_qagp(&func, breakpoints, nbreakpoints,
1e-8, 1e-8, n, &workspace, &result, &abserr);
n*=2;
if(ret!=GSL_SUCCESS) fprintf(stderr,"GSL error %s, increasing n to %li\n",gsl_strerror(ret),n);
gsl_set_error_handler(old_handler);
//gsl_set_error_handler(old_handler);
/* FIXME: do we care to keep the error estimate around? */
}
......@@ -269,14 +269,15 @@ log_radial_integrator *log_radial_integrator_init(double r1, double r2, int k, i
int interrupted=0;
OMP_BEGIN_INTERRUPTIBLE
integrator = malloc(sizeof(*integrator));
/* Temporarily turn off gsl_error handler which isn't thread safe. */
gsl_error_handler_t *old_handler = gsl_set_error_handler_off();
#pragma omp parallel for
for (size_t i = 0; i < size * size; i ++)
{
/*
if (OMP_WAS_INTERRUPTED)
OMP_EXIT_LOOP_EARLY;
*/
const size_t ix = i / size;
const size_t iy = i % size;
......@@ -288,7 +289,8 @@ log_radial_integrator *log_radial_integrator_init(double r1, double r2, int k, i
/* Note: using this where p > r0; could reduce evaluations by half */
z0[ix*size + iy] = log_radial_integral(r1, r2, p, b, k, cosmology);
}
gsl_set_error_handler(old_handler);
if (OMP_WAS_INTERRUPTED)
goto done;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment