Working up to normalisation

lalinference/src/LALInferenceInitCBC.c

@@ -1295,6 +1295,8 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
       /* If using margdist, remove the distance parameters and add the ranges into the model params as a way of passing them in */
       REAL8 a = log(Dmin), b=log(Dmax);
       LALInferenceAddMinMaxPrior(model->params, "logdistance", &a, &b, LALINFERENCE_REAL8_t);
+      UINT4 margdist=1;
+      LALInferenceAddVariable(model->params, "MARGDIST", &margdist, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);
 	  LALInferenceRemoveVariable(model->params, "logdistance");
   }
   

lalinference/src/LALInferenceLikelihood.c

@@ -47,11 +47,12 @@
 
 typedef enum
 {
-  GAUSSIAN,
-  STUDENTT,
-  MARGPHI,
-  MARGTIME,
-  MARGTIMEPHI
+  GAUSSIAN = 1,
+  STUDENTT = 2,
+  MARGPHI = 4,
+  MARGTIME = 8,
+  MARGTIMEPHI = 16,
+  MARGDIST = 32
 } LALInferenceLikelihoodFlags;
 
 
@@ -366,7 +367,14 @@ static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *cur
 
   /* ROQ likelihood stuff */
   REAL8 d_inner_h=0.0;
-
+  double dist_min, dist_max;
+  
+  UINT4 margdist = 0;
+  if(LALInferenceCheckVariable(model->params, "MARGDIST") && LALInferenceGetVariable(model->params, "MARGDIST"))
+  {
+      margdist = 1;
+      LALInferenceGetMinMaxPrior(model->params, "logdistance", &dist_min, &dist_max);
+  }
 
   if (LALInferenceCheckVariable(currentParams, "spcal_active") && (*(UINT4 *)LALInferenceGetVariable(currentParams, "spcal_active"))) {
     spcal_active = 1;
@@ -1121,6 +1129,12 @@ static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *cur
       /* This is marginalised over phase only for now */
       loglikelihood += -(S+D) + log(I0x) + R ;
       d_inner_h= 0.5*R;
+      
+      if(margdist )
+      {
+        loglikelihood = LALInferenceMarginalDistanceLogLikelihood(dist_min, dist_max, sqrt(S), d_inner_h);
+        loglikelihood -= D ;
+      }
       break;
     }
     case GAUSSIAN:
@@ -1163,8 +1177,16 @@ static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *cur
                   angMax = atan2(dh_S_phase->data[i], dh_S->data[i]);
                   xMax=x;
               }
-              double I0=log(gsl_sf_bessel_I0_scaled(x)) + fabs(x);
-              dh_S->data[i] = I0;
+              
+              if(margdist)
+              {
+                dh_S->data[i]=LALInferenceMarginalDistanceLogLikelihood(dist_min, dist_max, sqrt(S), x) - D;
+              }
+              else
+              {
+                double I0=log(gsl_sf_bessel_I0_scaled(x)) + fabs(x);
+                dh_S->data[i] = I0;
+              }
           }
       }
       size_t imax;
@@ -1208,16 +1230,27 @@ static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *cur
   if (OptimalSNR > 0.)
       MatchedFilterSNR = 2.0*d_inner_h/OptimalSNR;
     
-    if(1){
-        double dist_min, dist_max;
+
+    
+  LALInferenceAddVariable(currentParams,"optimal_snr",&OptimalSNR,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
+  LALInferenceAddVariable(currentParams,"matched_filter_snr",&MatchedFilterSNR,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
+
+  //loglikelihood = -1.0 * chisquared; // note (again): the log-likelihood is unnormalised!
+
+  return(loglikelihood);
+}
+
+double LALInferenceMarginalDistanceLogLikelihood(double dist_min, double dist_max, double OptimalSNR, double d_inner_h)
+{
 		static const size_t default_log_radial_integrator_size = 400;
-        LALInferenceGetMinMaxPrior(model->params, "logdistance", &dist_min, &dist_max);
-		static log_radial_integrator *integrator=NULL;
-		#pragma omp threadprivate(integrator)
+        double loglikelihood=0;
+		static log_radial_integrator *integrator;
+#pragma omp threadprivate(integrator)
+        double pmax = 100000; /* CHECKME: Max SNR allowed ? */
 
 		if (integrator == NULL)
 		{
-				double pmax = 20000; /* CHECKME: Max SNR allowed ? */
+                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(
@@ -1234,17 +1267,17 @@ static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *cur
 		
 
         //double marg_l = dist_integral(OptimalSNR*OptimalSNR, 2.0*d_inner_h, exp(dist_min), exp(dist_max));
-		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 = -D + marg_l;
-    }
-
-    
-  LALInferenceAddVariable(currentParams,"optimal_snr",&OptimalSNR,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
-  LALInferenceAddVariable(currentParams,"matched_filter_snr",&MatchedFilterSNR,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
-
-  //loglikelihood = -1.0 * chisquared; // note (again): the log-likelihood is unnormalised!
-
-  return(loglikelihood);
+        if (isnan(OptimalSNR) || isnan(d_inner_h) || pmax<OptimalSNR / sqrt(2))
+        {
+            loglikelihood=-INFINITY;
+            fprintf(stderr,"warning: Optimal SNR %lf exceeded pmax %lf\n",OptimalSNR/sqrt(2.0), pmax);
+        }
+        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;
+        }
+        return (loglikelihood);
 }
 
 /***************************************************************/

lalinference/src/LALInferenceLikelihood.h

@@ -179,6 +179,9 @@ REAL8 LALInferenceMarginalisedPhaseLogLikelihood(LALInferenceVariables *currentP
 
 REAL8 LALInferenceMarginalisedTimePhaseLogLikelihood(LALInferenceVariables *currentParams, LALInferenceIFOData *data, LALInferenceModel *model);
 
+/** Compute delta-log-likelihood for given distance min, max and OptimalSNR and d_inner_h when evaluated at 1Mpc */
+double LALInferenceMarginalDistanceLogLikelihood(double dist_min, double dist_max, double OptimalSNR, double d_inner_h);
+
 
 /**
  * Returns the log-likelihood marginalised over the time dimension

lalinference/src/distance_integrator.c

@@ -190,10 +190,10 @@ double log_radial_integral(double r1, double r2, double p, double b, int k, int
         log_offset = 0;
 
     params.scale = -log_offset;
-
-    {
+    size_t n = 64;
+    do {
         /* Maximum number of subdivisions for adaptive integration. */
-        static const size_t n = 64;
+
 
         /* Allocate workspace on stack. Hopefully, a little bit faster than
          * using the heap in multi-threaded code. */
@@ -216,13 +216,18 @@ 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();
 
         /* Perform adaptive Gaussian quadrature. */
         ret = gsl_integration_qagp(&func, breakpoints, nbreakpoints,
-            DBL_MIN, 1e-8, n, &workspace, &result, &abserr);
+            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);
 
         /* FIXME: do we care to keep the error estimate around? */
     }
+    while(ret!=GSL_SUCCESS);
 
     /* FIXME: do something with ret */
     (void)ret;
@@ -250,13 +255,15 @@ log_radial_integrator *log_radial_integrator_init(double r1, double r2, int k, i
     
     double *z1=calloc(size,sizeof(*z1));
     double *z2=calloc(size,sizeof(*z2));
-    double **z0=calloc(size,sizeof(*z0));
+    double *z0=calloc(size*size,sizeof(*z0));
     assert(z0 && z1 && z2);
+    /*
     for (size_t i=0;i<size;i++)
     {
         z0[i]=calloc(size,sizeof(*z0[i]));
         assert(z0[i]);
     }
+    */
     /* const double umax = xmax - vmax; */ /* unused */
 
     int interrupted=0;
@@ -279,27 +286,26 @@ log_radial_integrator *log_radial_integrator_init(double r1, double r2, int k, i
         const double r0 = exp(y);
         const double b = 2 * gsl_pow_2(p) / r0;
         /* Note: using this where p > r0; could reduce evaluations by half */
-        z0[ix][iy] = log_radial_integral(r1, r2, p, b, k, cosmology);
+        z0[ix*size + iy] = log_radial_integral(r1, r2, p, b, k, cosmology);
     }
 
 	if (OMP_WAS_INTERRUPTED)
         goto done;
 
-    region0 = bicubic_interp_init(*z0, size, size, xmin, ymin, d, d);
+    region0 = bicubic_interp_init(z0, size, size, xmin, ymin, d, d);
 
     for (size_t i = 0; i < size; i ++)
-        z1[i] = z0[i][size - 1];
+        z1[i] = z0[i*size + (size - 1)];
     region1 = cubic_interp_init(z1, size, xmin, d);
 
     for (size_t i = 0; i < size; i ++)
-        z2[i] = z0[i][size - 1 - i];
+        z2[i] = z0[i*size + (size - 1 - i)];
     region2 = cubic_interp_init(z2, size, umin, d);
 
 done:
     interrupted = OMP_WAS_INTERRUPTED;
     OMP_END_INTERRUPTIBLE
     
-    for (size_t i=0;i<size;i++) free(z0[i]);
     free(z2); free(z1); free(z0);
     if (interrupted || !(integrator && region0 && region1 && region2))
     {