add generation of distance posterior from marginalized likelihood

bilby/gw/conversion.py

 from __future__ import division
 import numpy as np
+from pandas import DataFrame
 
 from ..core.utils import logger, solar_mass
-from ..core.prior import DeltaFunction
+from ..core.prior import DeltaFunction, Interped
 
 try:
     from astropy.cosmology import z_at_value, Planck15
@@ -621,6 +622,11 @@ def _generate_all_cbc_parameters(sample, defaults, base_conversion,
     output_sample, _ = base_conversion(output_sample)
     output_sample = generate_mass_parameters(output_sample)
     output_sample = generate_spin_parameters(output_sample)
+    if likelihood is not None:
+        if likelihood.distance_marginalization:
+            output_sample = \
+                generate_distance_samples_from_marginalized_likelihood(
+                    output_sample, likelihood)
     output_sample = generate_source_frame_parameters(output_sample)
     compute_snrs(output_sample, likelihood)
     return output_sample
@@ -947,3 +953,42 @@ def compute_snrs(sample, likelihood):
 
     else:
         logger.debug('Not computing SNRs.')
+
+
+def generate_distance_samples_from_marginalized_likelihood(sample, likelihood):
+    if isinstance(sample, dict):
+        pass
+    elif isinstance(sample, DataFrame):
+        for ii in range(len(sample)):
+            temp = _generate_distance_sample_from_marginalized_likelihood(
+                dict(sample.iloc[ii]), likelihood)
+            sample['luminosity_distance'][ii] = temp['luminosity_distance']
+    return sample
+
+
+def _generate_distance_sample_from_marginalized_likelihood(sample, likelihood):
+    signal_polarizations = \
+        likelihood.waveform_generator.frequency_domain_strain(sample)
+    rho_mf_sq = 0
+    rho_opt_sq = 0
+    for ifo in likelihood.interferometers:
+        signal = ifo.get_detector_response(signal_polarizations, sample)
+        rho_mf_sq += ifo.matched_filter_snr_squared(signal=signal)
+        rho_opt_sq += ifo.optimal_snr_squared(signal=signal)
+
+    rho_mf_sq_dist = \
+        rho_mf_sq * sample['luminosity_distance'] / \
+        likelihood._distance_array
+
+    rho_opt_sq_dist = \
+        rho_opt_sq * sample['luminosity_distance']**2 / \
+        likelihood._distance_array**2
+
+    distance_log_like = (rho_mf_sq_dist.real - rho_opt_sq_dist.real / 2)
+
+    distance_post = np.exp(distance_log_like - max(distance_log_like)) *\
+        likelihood.distance_prior_array
+
+    sample['luminosity_distance'] = Interped(
+        likelihood._distance_array, distance_post).sample()
+    return sample