Calibration plotting

bilby/core/utils.py

@@ -703,6 +703,31 @@ def logtrapzexp(lnf, dx):
     return np.log(dx / 2.) + logsumexp([logsumexp(lnf[:-1]), logsumexp(lnf[1:])])
 
 
+def credible_interval(samples, confidence_level=.9, lower=True):
+    """
+    Return location of lower or upper confidence levels
+    Based on lalinference.bayespputils.cred_interval
+
+    Parameters
+    ----------
+    x: List of samples.
+    cl: Confidence level to return the bound of.
+    lower: If ``True``, return the lower bound, otherwise return the upper bound.
+
+    Returns
+    -------
+    float: the upper or lower confidence level
+
+    """
+    def cred_level(cl, x):
+        return np.sort(x, axis=0)[int(cl * len(x))]
+
+    if lower:
+        return cred_level((1. - confidence_level) / 2, samples)
+    else:
+        return cred_level((1. + confidence_level) / 2, samples)
+
+
 def run_commandline(cl, log_level=20, raise_error=True, return_output=True):
     """Run a string cmd as a subprocess, check for errors and return output.
 

bilby/gw/result.py

 from __future__ import division
 
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+
 from ..core.result import Result as CoreResult
 from ..core.utils import logger
+from .utils import (plot_spline_pos,
+                    spline_angle_xform)
 
 
 class CompactBinaryCoalesenceResult(CoreResult):
@@ -80,5 +86,80 @@ class CompactBinaryCoalesenceResult(CoreResult):
             logger.info("No injection for detector {}".format(detector))
             return None
 
+    def plot_calibration_posterior(self, level=.9):
+        """ Plots the calibration amplitude and phase uncertainty.
+        Adapted from the LALInference version in bayespputils
+
+        Parameters
+        ----------
+        level: float,  percentage for confidence levels
+
+        Returns
+        -------
+        saves a plot to outdir+label+calibration.png
+
+        """
+        fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(15, 15), dpi=500)
+        posterior = self.posterior
+
+        font_size = 32
+        outdir = self.outdir
+
+        parameters = posterior.keys()
+        ifos = np.unique([param.split('_')[1] for param in parameters if 'recalib_' in param])
+        if ifos.size == 0:
+            logger.info("No calibration parameters found. Aborting calibration plot.")
+            return
+
+        for ifo in ifos:
+            if ifo == 'H1':
+                color = 'r'
+            elif ifo == 'L1':
+                color = 'g'
+            elif ifo == 'V1':
+                color = 'm'
+            else:
+                color = 'c'
+
+            # Assume spline control frequencies are constant
+            freq_params = np.sort([param for param in parameters if
+                                   'recalib_{0}_frequency_'.format(ifo) in param])
+
+            logfreqs = np.log([posterior[param].iloc[0] for param in freq_params])
+
+            # Amplitude calibration model
+            plt.sca(ax1)
+            amp_params = np.sort([param for param in parameters if
+                                  'recalib_{0}_amplitude_'.format(ifo) in param])
+            if len(amp_params) > 0:
+                amplitude = 100 * np.column_stack([posterior[param] for param in amp_params])
+                plot_spline_pos(logfreqs, amplitude, color=color, level=level,
+                                label="{0} (mean, {1}%)".format(ifo.upper(), int(level * 100)))
+
+            # Phase calibration model
+            plt.sca(ax2)
+            phase_params = np.sort([param for param in parameters if
+                                    'recalib_{0}_phase_'.format(ifo) in param])
+            if len(phase_params) > 0:
+                phase = np.column_stack([posterior[param] for param in phase_params])
+                plot_spline_pos(logfreqs, phase, color=color, level=level,
+                                label="{0} (mean, {1}%)".format(ifo.upper(), int(level * 100)),
+                                xform=spline_angle_xform)
+
+        ax1.tick_params(labelsize=.75 * font_size)
+        ax2.tick_params(labelsize=.75 * font_size)
+        plt.legend(loc='upper right', prop={'size': .75 * font_size}, framealpha=0.1)
+        ax1.set_xscale('log')
+        ax2.set_xscale('log')
+
+        ax2.set_xlabel('Frequency (Hz)', fontsize=font_size)
+        ax1.set_ylabel('Amplitude (%)', fontsize=font_size)
+        ax2.set_ylabel('Phase (deg)', fontsize=font_size)
+
+        filename = os.path.join(outdir, self.label + '_calibration.png')
+        fig.tight_layout()
+        fig.savefig(filename, bbox_inches='tight')
+        plt.close(fig)
+
 
 CBCResult = CompactBinaryCoalesenceResult

bilby/gw/utils.py

@@ -3,10 +3,13 @@ import os
 import json
 
 import numpy as np
+from scipy import interpolate
+import matplotlib.pyplot as plt
 
 from ..core.utils import (gps_time_to_gmst, ra_dec_to_theta_phi,
                           speed_of_light, logger, run_commandline,
-                          check_directory_exists_and_if_not_mkdir)
+                          check_directory_exists_and_if_not_mkdir,
+                          credible_interval)
 
 try:
     from gwpy.timeseries import TimeSeries
@@ -856,3 +859,69 @@ def lalsim_SimInspiralWaveformParamsInsertTidalLambda2(
 
     return lalsim.SimInspiralWaveformParamsInsertTidalLambda2(
         waveform_dictionary, lambda_2)
+
+
+def spline_angle_xform(delta_psi):
+    """
+    Returns the angle in degrees corresponding to the spline
+    calibration parameters delta_psi.
+    Based on the same function in lalinference.bayespputils
+
+    Parameters
+    ----------
+    delta_psi: calibration phase uncertainity
+
+    Returns
+    -------
+    float: delta_psi in degrees
+
+    """
+    rotation = (2.0 + 1.0j * delta_psi) / (2.0 - 1.0j * delta_psi)
+
+    return 180.0 / np.pi * np.arctan2(np.imag(rotation), np.real(rotation))
+
+
+def plot_spline_pos(log_freqs, samples, nfreqs=100, level=0.9, color='k', label=None, xform=None):
+    """
+    Plot calibration posterior estimates for a spline model in log space.
+    Adapted from the same function in lalinference.bayespputils
+
+    Parameters
+    ----------
+        log_freqs: The (log) location of spline control points.
+        samples: List of posterior draws of function at control points ``log_freqs``
+        nfreqs: Number of points to evaluate spline at for plotting.
+        level: Credible level to fill in.
+        color: Color to plot with.
+        label: Label for plot.
+        xform: Function to transform the spline into plotted values.
+
+    """
+    freq_points = np.exp(log_freqs)
+    freqs = np.logspace(min(log_freqs), max(log_freqs), nfreqs, base=np.exp(1))
+
+    data = np.zeros((samples.shape[0], nfreqs))
+
+    if xform is None:
+        scaled_samples = samples
+    else:
+        scaled_samples = xform(samples)
+
+    mu = np.mean(scaled_samples, axis=0)
+    lower_confidence_level = mu - credible_interval(scaled_samples, level, lower=True)
+    upper_confidence_level = credible_interval(scaled_samples, level, lower=False) - mu
+    plt.errorbar(freq_points, mu, yerr=[lower_confidence_level, upper_confidence_level],
+                 fmt='.', color=color, lw=4, alpha=0.5, capsize=0)
+
+    for i, sample in enumerate(samples):
+        temp = interpolate.spline(log_freqs, sample, np.log(freqs))
+        if xform is None:
+            data[i] = temp
+        else:
+            data[i] = xform(temp)
+
+    line, = plt.plot(freqs, np.mean(data, axis=0), color=color, label=label)
+    color = line.get_color()
+    plt.fill_between(freqs, credible_interval(data, level), credible_interval(data, level, lower=False),
+                     color=color, alpha=.1, linewidth=0.1)
+    plt.xlim(freq_points.min() - .5, freq_points.max() + 50)